In [44]:
from IPython.display import HTML

HTML('''<script>
code_show=true; 
function code_toggle() {
 if (code_show){
 $('div.input').hide();
 } else {
 $('div.input').show();
 }
 code_show = !code_show
} 
$( document ).ready(code_toggle);
</script>
<form action="javascript:code_toggle()"><input type="submit" value="Click here to toggle on/off the raw code."></form>''')
Out[44]:
In [2]:
%%html
<style>
td {
  font-size: 14px
}
table {
  float:left
}
</style>

NSBH Spin Analysis


In this work, I consider mass ratio 5 NSBH mergers with varying black hole spin. Neutron star is non-spinning with compactness of 0.15. Black Hole spin varies from -0.5 to 0.75. As a test study, the initial separation is kept small (only 7M). Below table provides additional simulation details.

Simulation Type Mass Ratio Separation BH Spin Resolution TOV Resolution BH1
Sim1 NSBH 5 7M -0.5 $R_*$/51 $M_{BH}$/90
Sim2 NSBH 5 7M 0.0 $R_*$/50 $M_{BH}$/89
Sim3 NSBH 5 7M 0.25 $R_*$/50 $M_{BH}$/89
Sim4 NSBH 5 7M 0.5 $R_*$/50 $M_{BH}$/90
Sim5 NSBH 5 7M 0.75 $R_*$/50 $M_{BH}$/89
In [3]:
#Import Statements
import numpy as np
import matplotlib.pyplot as plt
import os, glob
from itertools import cycle
from shutil import copyfile
from scipy.interpolate import splrep, splev


import yt
from yt.visualization.eps_writer import multiplot_yt

from mpl_toolkits.axes_grid1 import AxesGrid
yt.funcs.mylog.setLevel(50)


import matplotlib as mpl
mpl.rc('lines', linewidth=2, color='r')
mpl.rc('font', size=18)
mpl.rc('axes', labelsize=18, grid=True)
mpl.rc('xtick', labelsize=14)
mpl.rc('ytick', labelsize=14)
mpl.rc('legend', fontsize=16)
/localdata/bkhamesra3/softwares/anaconda3/envs/py27/lib/python2.7/site-packages/h5py/__init__.py:36: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.
  from ._conv import register_converters as _register_converters
In [4]:
from scipy.constants import G, c
Msun = 1.989e30
m1_in_Msun = c**2/G/Msun
s1_in_Msun = c**3/G/Msun

#print("1m = %g Msun, 1s = %g Msun" %(m1_in_Msun, s1_in_Msun))
In [5]:
#Make line plots, log and semilog plots of one variable
def lineplots(variable_name, filename, grid_coord='x',xlimit=None, ylimit=None, shiftx=0.,shifty=0., logx=False, logy=False, num_snapshots=None, title=None ):
    ''' Function - Make line plots, log and semilog plots of one variable
        variable - string input - variable name
        filename - string input - filepath
        xlimit - optional - float - limit of x axis
        logx - binary (True/False) - to change scale of x axis to log
        logy - binary (True/False) - to change scale of y axis to log
        title - title of the plot
        num_snapshots - number of snapshots in time of variable
        shiftx - shift x coordinate by this amount
        
        
    '''
    if grid_coord=='x': grid_coln = 9
    elif grid_coord=='y': grid_coln = 10
    else: grid_coln=11
    
    t, x, var = np.loadtxt(filename, unpack=True, usecols=(8,grid_coln,12))
    
    t_uniq = np.unique(t)
   
    plt.figure(figsize=(15,8))
    if num_snapshots==None:
        for time_instance in t_uniq:
            
            xx_instance = x[t==time_instance]
            var_instance = var[t==time_instance]
            
            idx_sorted = np.argsort(xx_instance)
            xx_sort, var_sort = xx_instance[idx_sorted], var_instance[idx_sorted]
            
            legend=False
            plt.plot(xx_sort,var_sort)
      
    else:
        legend=True
        time_snapshot_idx = np.linspace(0,len(t_uniq)-2, num_snapshots)
        time_snapshot_idx = time_snapshot_idx.astype(int)
     
        time_snapshot = t_uniq[time_snapshot_idx]
        
        for time_instance in time_snapshot:
            idx_instance =  np.where(t==time_instance)
            
            xx_instance = x[t==time_instance]
            var_instance = var[t==time_instance]
            
            idx_sorted = np.argsort(xx_instance)
            xx_sort, var_sort = xx_instance[idx_sorted], var_instance[idx_sorted]
            plt.plot(xx_sort,var_sort, label="t=%g"%time_instance)
        
    plt.xlabel('x')
    plt.ylabel(variable_name)
    if not(title==None):
        plt.title(title, fontsize=16)
    if isinstance(xlimit,float):
        plt.xlim(-1*xlimit+shiftx,xlimit+shiftx )
    if isinstance(ylimit,float):
        plt.ylim(-1*ylimit+shifty,ylimit+shifty )
    
    if(logx==True): plt.xscale('log')
    if(logy==True): plt.yscale('log')
    if legend==True:
        plt.legend()
    plt.show()
    plt.close()
    



def linesubplots(numvars, filename, yaxisname, labelname):
    ''' Function - Make subplots of multiple variables
    numvars - integer - number of variables 
    filename - list of file paths
    yaxisname - list of y-axis labels
    labelname - list of labels 
    '''
    
    num_rows, num_cols = 0, 0
    if numvars%2==0:
        num_rows, num_cols = numvars/2, 2
    else:
        num_rows, num_cols = numvars, 1
        
    fig = plt.figure(figsize=(12,7*num_rows))
    
    time, x, variable = {}, {}, {}
    
    #Read data and store in dictionaries
    for i in range(0,numvars):
        tt, xx, data = np.loadtxt(filename[i], unpack=True, usecols=(8,9,12))
        idx_sorted = np.argsort(xx)
        xx, data = xx[idx_sorted], data[idx_sorted]
    
    
        time['var%d'%i] = tt
        x['var%d'%i] = xx
        variable['var%d'%i] = data
     
        plt.subplot(num_rows, num_cols, i+1)
        plt.plot(x['var%d'%i], variable['var%d'%i], label=labelname[i])
        plt.ylabel(yaxisname[i])
        plt.xlabel('x')
        plt.xlim(-8,8)
        plt.legend()
    
    plt.tight_layout()    
    plt.show()
    plt.close()
        
    
    
        
def comparisonplots(varname, numdatasets, filename, labelname, xlimit=None, logx=False, logy=False, freq=1):
    ''' Function - single plots showing comparison of two or more variables 
    varname - string - variable name which is being compared
    numdatasets - integer - number of datasets 
    filename - list of file paths 
    labelname - list of label of each dataset
    '''
    time, x, variable = {}, {}, {}
    lines = ["-","--","-.",":"]
    linecycler = cycle(lines)
    fig = plt.figure(figsize=(15,8))
   
    #Read data and store in dictionaries
    for i in range(numdatasets):
        tt, xx, data = np.loadtxt(filename[i], unpack=True, usecols=(8,9,12))
        t_arr = np.unique(tt)

        linestyle=next(linecycler)
        for t_elem in t_arr: 
            xx_t = xx[np.where(tt==t_elem)]
            data_t = data[np.where(tt==t_elem)]
            idx_sorted = np.argsort(xx_t)
            xx_t, data_t = xx_t[idx_sorted], data_t[idx_sorted]
           
            plt.plot(xx_t, data_t, ls=linestyle, label=labelname[i])
    
    plt.ylabel(varname)
    plt.xlabel('x')
    if(logx==True): plt.xscale('log')
    if(logy==True): plt.yscale('log')
   
    if isinstance(xlimit,float):
        plt.xlim(-1*xlimit,xlimit)
    plt.legend()
  
    plt.show()
    plt.close()
        
    

def comparisonmultiplots(varnames, filenames, labelnames, xlimit=None, logx=False, logy=False):
    ''' Create multiple subplots comparing two quantities in each subplot
    varnames: list of name of each variable (this name is also used to label y axis)
    numdatasets: number of datasets for each variable
    filename: dictionary of filenames for each variable and each dataset. Use index in form of rows and columns where 
    each row contains all datasets for one subplot
    labelname: list of names of labels for each varname. Should have size of numdatasets
    '''

    num_vars = len(varnames)
    num_datasets = len(labelnames)
    num_files = len(filenames.keys())
    assert (num_files==num_vars*num_datasets), "Total number of filenames do not match with number of variables and labelnames"
    
    if xlimit==None:
        xlimit=30
        
    num_cols=2
    if num_vars%2==0:
        num_rows = num_vars/2    
    else:
        num_rows, num_cols = num_vars, 1
        
    
    
    time, x, variable={},{},{}
    plt.figure(figsize=(15, 8*num_rows))
    linestyles = ["-","--","-.",":"]
    
    for i in range(num_files):
        
        tt, xx, data = np.loadtxt(filenames[sorted(filenames.keys())[i]], unpack=True, usecols=(8,9,12))
        idx_sorted = np.argsort(xx)
        xx, data = xx[idx_sorted], data[idx_sorted]
    
        time[i] = tt
        x[i] = xx
        variable[i] = data
        
        if i%num_datasets==0:
            
            plt.subplot(num_rows, num_cols, 1+i/num_datasets)
            plt.ylabel(varnames[i/num_datasets])
            plt.xlabel('x')
            plt.xlim(-1*xlimit, xlimit)
            
        
        #print i
        plt.plot(x[i], variable[i], ls = linestyles[i%num_datasets], label=labelnames[i%num_datasets])
        plt.legend(fontsize=14)
    plt.show()
    plt.close()
        
        
In [6]:
os.chdir('/localdata2/bkhamesra3/simulations/Hive/NSBH/MassRatio_5/SpinStudy/')
In [7]:
#Directories
nsbh_q5_d7_a1_n0d5_dir1  = 'NSBH_q5_D7_a1_-0.5_a2_0.0_C0.15_M109/'
nsbh_q5_d7_a1_0d0_dir2   = 'NSBH_q5_D7_a1_0.0_a2_0.0_C0.15_M109/'
nsbh_q5_d7_a1_p0d25_dir3 = 'NSBH_q5_D7_a1_0.25_a2_0.0_C0.15_M107/'
nsbh_q5_d7_a1_p0d5_dir4  = 'NSBH_q5_D7_a1_0.5_a2_0.0_C0.15_M109/'
nsbh_q5_d7_a1_p0d75_dir5 = 'NSBH_q5_D7_a1_0.75_a2_0.0_C0.15_M107/'


Madm_nsbh_a1_0d0 =   0.9938816335979512306
Madm_nsbh_a1_p0d25 = 0.9937419751277200008
Madm_nsbh_a1_p0d5 =  0.993621911714029471
Madm_nsbh_a1_p0d75 = 0.9935297248578904838
Madm_nsbh_a1_n0d5 =  0.9942234586460986234



#Initial Configuration
Dist_from_src=100. #Mpc
Mtotal = 1.35*3  #Msun
mass_ratio=5
MPc = 3.086e+22  #metres
Mtotal_to_MPc = (Mtotal/m1_in_Msun/MPc)

nsbh_a1_n0d5 = {'dir':nsbh_q5_d7_a1_n0d5_dir1, 'label':r'$a_{BH} = -0.5$', 'type':'NSBH', 'symmetry':True, \
                'ls':'-', 'color':'k', 'M_adm':Madm_nsbh_a1_n0d5} #,  'tidal_radius':6.6} 

nsbh_a1_0d0 = {'dir':nsbh_q5_d7_a1_0d0_dir2 , 'label':r'$a_{BH} = 0$', 'type':'NSBH', 'symmetry':True, \
               'ls':':', 'color':'b', 'M_adm':Madm_nsbh_a1_0d0}

nsbh_a1_0d25 = {'dir':nsbh_q5_d7_a1_p0d25_dir3, 'label':r'$a_{BH} = 0.25$' ,'type':'NSBH', 'symmetry':True, \
               'ls':':', 'color':'orangered', 'M_adm':Madm_nsbh_a1_p0d25}

nsbh_a1_0d5 = {'dir':nsbh_q5_d7_a1_p0d5_dir4, 'label':r'$a_{BH} = 0.5$', 'type':'NSBH', 'symmetry':True, \
               'ls':':', 'color':'m', 'M_adm':Madm_nsbh_a1_p0d5}

nsbh_a1_0d75 = {'dir':nsbh_q5_d7_a1_p0d75_dir5, 'label':r'$a_{BH} = 0.75$', 'type':'NSBH', 'symmetry':True, \
               'ls':':', 'color':'forestgreen', 'M_adm':Madm_nsbh_a1_p0d75}


compare_nsbh = {'nsbh1':nsbh_a1_n0d5, 'nsbh2':nsbh_a1_0d0, 'nsbh3':nsbh_a1_0d25, 'nsbh4':nsbh_a1_0d5, 'nsbh5':nsbh_a1_0d75}
figdir = '/localdata2/bkhamesra3/simulations/Hive/NSBH/Spin_Comparison/Sep_7M'

if not os.path.exists(figdir):
    os.makedirs(figdir)

Merger Time


  • $T_{Peak}$ - Time of $\psi_4$ peak amplitude.
  • $T_{Peak} - R_{ex}$ - Approximate Retarded Time of peak amplitude defined as $t' = t - R_{ex}$. This is used as the definition of the Time of Merger in this analysis.
  • $T_{Peak}$ (Retarded) - Retarded Time of peak amplitude defined as - $t' = t - R_{ex} - 2*log(\frac{R_{ex}}{2}-1)$
In [8]:
# Define Merger Time in retarded frame - t' = t - D - 2*log(D/2-1) - 
# Reference - Equation 26, https://iopscience.iop.org/article/10.1088/0264-9381/31/2/025012/pdf

print("Simulation \t \tT_peak \t T_peak - R_ex \t T_peak(retarded)")
print("-"*100)
for i, key in enumerate(sorted(compare_nsbh.keys())):

    if (compare_nsbh[key]['type']=='NSBH'):
        datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
    elif(compare_nsbh[key]['type']=='BBH'):
        datadir = compare_nsbh[key]['dir']

    psi4 = os.path.join(datadir, 'Ylm_WEYLSCAL4::Psi4r_l2_m2_r70.00.asc')
    if not os.path.exists(psi4):
        psi4 = os.path.join(datadir, 'mp_Psi4r_l2_m2_r70.00.asc')
    
    t, re, im = np.loadtxt(psi4, unpack=True, usecols=(0,1,2))
    det_radius = float(((psi4.split('_')[-1]).split('.asc')[0]).split('r')[1])
    amp = np.sqrt(re**2 + im**2)
    tmax_det = t[amp==np.max(amp)]
    tmax_ret = tmax_det-det_radius -2.*np.log(det_radius/2. - 1.)
    tmerger_src = tmax_det-det_radius 
    compare_nsbh[key]['t_retarded'] = tmax_ret
    compare_nsbh[key]['t_merger_src'] = tmerger_src
    print("%15s \t %.2f \t %.2f \t %.2f" %(compare_nsbh[key]['label'],tmax_det, compare_nsbh[key]['t_merger_src'],compare_nsbh[key]['t_retarded']))
    
Simulation 	 	T_peak 	 T_peak - R_ex 	 T_peak(retarded)
----------------------------------------------------------------------------------------------------
$a_{BH} = -0.5$ 	 225.62 	 155.62 	 148.56
   $a_{BH} = 0$ 	 399.32 	 329.32 	 322.26
$a_{BH} = 0.25$ 	 523.26 	 453.26 	 446.21
 $a_{BH} = 0.5$ 	 653.93 	 583.93 	 576.88
$a_{BH} = 0.75$ 	 784.30 	 714.30 	 707.25
In [9]:
#Compute ADM angular momentum
for i, key in enumerate(sorted(compare_nsbh.keys())):

    if (compare_nsbh[key]['type']=='NSBH'):
        datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
    
    parameter_file = glob.glob(os.path.join(datadir, '*.par'))[0]
    labelname = compare_nsbh[key]['label']
    #adm_angmom = compare_nsbh[key]['J_adm']
    
    
    x1,y1,z1,x2,y2,z2,px1,py1,pz1,px2,py2,pz2,jx1,jy1,jz1,jx2,jy2,jz2 = np.zeros(18)
    
    
    with open(parameter_file) as parfile:
        for line in parfile.readlines():
            if  'BowenID::object_rx[0]' in line:
                x1 = float(line.split('=')[1])
            elif 'BowenID::object_rx[1]' in line:
                x2 = float(line.split('=')[1])
            elif 'BowenID::object_Px[0]' in line:
                px1 = float(line.split('=')[1])
            elif 'BowenID::object_Py[0]' in line:
                py1 = float(line.split('=')[1])
            elif 'BowenID::object_Pz[0]' in line:
                pz1 = float(line.split('=')[1])
            elif 'BowenID::object_Px[1]' in line:
                px2 = float(line.split('=')[1])
            elif 'BowenID::object_Py[1]' in line:
                py2 = float(line.split('=')[1])
            elif 'BowenID::object_Pz[1]' in line:
                pz2 = float(line.split('=')[1])
            elif 'BowenID::object_Jx[0]' in line:
                jx1 = float(line.split('=')[1])
            elif 'BowenID::object_Jy[0]' in line:
                jy1 = float(line.split('=')[1])
            elif 'BowenID::object_Jz[0]' in line:
                jz1 = float(line.split('=')[1])
            elif 'BowenID::object_Jx[1]' in line:
                jx2 = float(line.split('=')[1])
            elif 'BowenID::object_Jy[1]' in line:
                jy2 = float(line.split('=')[1])
            elif 'BowenID::object_Jz[1]' in line:
                jz2 = float(line.split('=')[1])
                
    r_bh = np.array((x1, y1, z1))
    p_bh = np.array((px1, py1, pz1))
    j_bh = np.array((jx1, jy1, jz1))
   
    r_ns = np.array((x2, y2, z2))
    p_ns = np.array((px2, py2, pz2))
    j_ns = np.array((jx2, jy2, jz2))
            
    j_adm = np.cross(r_bh, p_bh ) + np.cross(r_ns, p_ns) + j_bh + j_ns
    
    compare_nsbh[key]['J_adm'] = j_adm
    

Density Snapshots


We look at the density evolution of the NSBH systems. Figures from top to bottom are arranged in the increasing order of the BH spins. For each case, the panels on the top from left to right show time snapshots at 100M, 50M and 20M before the merger while panels on the bottom show the snapshots at the merger, and 100M and 300 M after the merger. A clear impact of the merger is the increasing disruption of the star which is minimal for $a_{BH} = -0.5$ and maximum for $a_{BH} = 0.75$.

In [10]:
#Define Fields for YT
from yt import derived_field 

def _density_gu(field, data):
    '''Specifies density in geometrical units'''
    return  data["density"]

def _density_cgs_q5(field, data):
    '''Specifies density in geometrical units'''
    Mt = 1.35*6
    return  (data["density"]/Mt**2)


#yt.add_field(("gas", "density_gu"), function=_density_gu, units="g/cm**3")
yt.add_field(("gas", "density_gu"), function=_density_gu, units="code_mass/code_length**3")
#yt.add_field(("gas", "density_cgs_q2"), function=_density_cgs_q2, units="g/cm**3")
#yt.add_field(("gas", "density_cgs_q3"), function=_density_cgs_q3, units="g/cm**3")
yt.add_field(("gas", "density_cgs_q5"), function=_density_cgs_q5, units="g/cm**3")


for i, key in enumerate(sorted(compare_nsbh.keys())):

    labelname = compare_nsbh[key]['label']
    #print(labelname+'\n')
    if (compare_nsbh[key]['type']=='NSBH'):
        datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
        t_merger_src = compare_nsbh[key]['t_merger_src']

       # print("%s:Time of merger = %g"%(labelname, t_merger_src))
    else:
        continue
        
    M_in_Msun  = 1.35*6
    fieldname = "density_cgs_q5"
        
    #filename = glob.glob(os.path.join(datadir,'../../output-0000/NSBH*/rho.x.asc'))[0]
    
    #t, x, rho = np.loadtxt(filename, unpack=True, usecols=(8,9,12))
    #x = x[t==0]
    #rho = rho[t==0]
    max_rho = 0.101473916805118 #np.amax(rho)
    #x_maxrho = x[np.argmax(rho)]
    max_rho_si = (max_rho/(M_in_Msun**2))*(Msun*m1_in_Msun**3)
    
    h5file = os.path.join(datadir, 'rho.xy.h5') 
    varname1 = fieldname
    ds1 = yt.load(h5file, iteration=0) # load data
    _, max_rho_si_yt = ds1.all_data().quantities.extrema(varname1)
    
    varname2 = 'density_gu'
    ds2 = yt.load(h5file, iteration=0) # load data
    min_rho_gu_yt, max_rho_gu_yt = ds2.all_data().quantities.extrema(varname2)

    #print("t=0: x=%g, rho_c = %g/M^2 (%g g/cm^3), rho_c (yt) = %g /M^2 (%g g/cm^3)"%(x_maxrho, max_rho, max_rho_si/1e3, max_rho_gu_yt, max_rho_si_yt))
        
Because 'sampling_type' not specified, yt will assume a cell 'sampling_type'
In [11]:
#Rest Mass Density Contour Plots

import yt
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import AxesGrid
mpl.rc('axes', labelsize=16, grid=False)


    # See http://matplotlib.org/mpl_toolkits/axes_grid/api/axes_grid_api.html
    # These choices of keyword arguments produce a four panel plot that includes
    # four narrow colorbars, one for each plot.  Axes labels are only drawn on the
    # bottom left hand plot to avoid repeating information and make the plot less
    # cluttered.
    
for i, key in enumerate(sorted(compare_nsbh.keys())):
    
    labelname = compare_nsbh[key]['label']
    spin = float(labelname.split('=')[1][:-1])
    
    if (compare_nsbh[key]['type']=='NSBH'):
        datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
        t_merger_src = compare_nsbh[key]['t_merger_src']
    else:
        continue
   
    fig = plt.figure(figsize=(10,8))

    #Setup axes grid
    grid = AxesGrid(fig, (0.1,0.1,1.2,1.2),
                    nrows_ncols = (2, 3),
                    axes_pad = 0.5,
                    label_mode = "L",
                    share_all = False,
                    cbar_location="right",
                    cbar_mode="single",
                    cbar_size="3%",
                    cbar_pad="0%")
    
    #Define variable name for plots
    varname = 'density'
    
    #Path for density and BH diagnostics data
    h5file = os.path.join(datadir, 'rho.xy.h5') 
    bhdiag = os.path.join(datadir, 'BH_diagnostics.ah1.gp')
    shifttracker = os.path.join(datadir, 'ShiftTracker0.asc')

    assert(os.path.exists(h5file)),"H5 file not found"


    M = 1.35*6
    #varname_cgs = "density_cgs_q5"
    
    #Read data from BH diagnostics file
    it_st, t_st = np.loadtxt(shifttracker, unpack=True, usecols=(0,1))
    it, t, xbh, ybh, zbh, rbh, bh_mass = np.loadtxt(bhdiag, unpack=True, usecols=(0,1,2,3,4,7,26))
    idx_merger = np.amin(np.where(t>=t_merger_src))
    
    #print(t_merger_src, t[idx_merger])
    t_centered = t - t[idx_merger]
    t_st_centered = t_st - t[idx_merger]
    t_premerger1 = -100
    t_premerger2=-50
    t_premerger3=-20
    t_postmerger1=0
    t_postmerger2=100
    t_postmerger3=300
    
    
    time_inst = [t_premerger1, t_premerger2, t_premerger3, t_postmerger1, t_postmerger2, t_postmerger3]
    ds_allit = yt.load(h5file)
    iterations = ds_allit.get_all_iterations()
    k=0
    
    print("BH Spin = {} \n\t Pre Merger - t1 = {}({}), t2 = {}({}), t3 = {}({}),\
          \n\t Post Merger - t1 = {}({}), t2 = {}({}), t3 = {}({})\n"\
          .format(spin, t_premerger1,t_premerger1 + t[idx_merger], t_premerger2, t_premerger2 + t[idx_merger],\
                   t_premerger3, t_premerger3 + t[idx_merger], t_postmerger1, t_postmerger1 + t[idx_merger],t_postmerger2,\
                  t_postmerger2 + t[idx_merger],t_postmerger3, t_postmerger3 + t[idx_merger]))
    
    for t_pm in time_inst:
        it_pm_st = it_st[np.amax(np.where(t_st_centered<=t_pm))]
        idx_pm = np.amax(np.where(t_centered<=t_pm))
        it_pm  = it [idx_pm]
        xbh_pm = xbh[idx_pm]
        ybh_pm = ybh[idx_pm]
        zbh_pm = zbh[idx_pm]
        rbh_pm = rbh[idx_pm]
        
        
        iter_premerger = np.amax(np.where(iterations<=it_pm))
        ds = yt.load(h5file, iteration=iterations[iter_premerger]) # load data
        
        min_var, max_var = ds.all_data().quantities.extrema(varname)
        
        #As we cannot change units while loading the data, we call for new variables with specific units
        #min_rho_gu, max_rho_gu = ds.all_data().quantities.extrema('density_gu')
        #min_rho_si, max_rho_si = ds.all_data().quantities.extrema(varname_cgs)
        #print("t=%d: Maximum Density = %g/M^2 (%g g/cm^3)"%(t_pm, max_rho_gu, max_rho_si))
        
        
        axis_width=30
        #Plotting density is fine here as we later set the units to code units which plots exact data stored in the file
        p = yt.SlicePlot(ds, 'z', varname, axes_unit='code_length', origin='native', width=axis_width, center = [0,0,0])
        p.set_unit('density', 'code_mass/code_length**3')
        p.set_font_size(16)
        p.set_xlabel('x/M')
        p.set_ylabel('y/M')
        
        p.set_colorbar_label('density','Density (1/M$^2$)')

        p.set_cmap(varname, 'dusk') #'viridis')
        p.annotate_sphere([xbh_pm, ybh_pm, zbh_pm], radius=(rbh_pm, 'code_length'), circle_args={'facecolor':'black','fill':1, 'edgecolor':'white'})
        #p.annotate_contour(varname, plot_args={'colors':'white', 'linestyles':'dotted'})
        p.set_font_size(20)
        #p.annotate_text((0.45,0.9), 't = %.2f'%t[idx_pm], coord_system='figure', text_args={'color':'white'})
        
        #if k==i:p.annotate_title('q=%d'%mass_ratio)
       
        p.set_zlim(varname, 0.05, 1e-7)
               
        
        plot = p.plots[varname]
                  
        plot.figure = fig
        plot.axes = grid[k].axes
        plot.cax = grid.cbar_axes[k]
    
        # Actually redraws the plot.
        p._setup_plots()
        k+=1
        
  
        #plt.savefig( os.path.join(figdir, 'DensitySnapshots_q%d.png')%mass_ratio, dpi=100 )
    plt.show()
The is_string_like function was deprecated in version 2.1.
The is_string_like function was deprecated in version 2.1.
BH Spin = -0.5 
	 Pre Merger - t1 = -100(55.616), t2 = -50(105.616), t3 = -20(135.616),          
	 Post Merger - t1 = 0(155.616), t2 = 100(255.616), t3 = 300(455.616)

BH Spin = 0.0 
	 Pre Merger - t1 = -100(229.327), t2 = -50(279.327), t3 = -20(309.327),          
	 Post Merger - t1 = 0(329.327), t2 = 100(429.327), t3 = 300(629.327)

BH Spin = 0.25 
	 Pre Merger - t1 = -100(353.264), t2 = -50(403.264), t3 = -20(433.264),          
	 Post Merger - t1 = 0(453.264), t2 = 100(553.264), t3 = 300(753.264)

BH Spin = 0.5 
	 Pre Merger - t1 = -100(483.942), t2 = -50(533.942), t3 = -20(563.942),          
	 Post Merger - t1 = 0(583.942), t2 = 100(683.942), t3 = 300(883.942)

BH Spin = 0.75 
	 Pre Merger - t1 = -100(614.299), t2 = -50(664.299), t3 = -20(694.299),          
	 Post Merger - t1 = 0(714.299), t2 = 100(814.299), t3 = 300(1014.299)

Density Oscillations__


Below figure shows the density oscillations in the neutron star and corresponding frequency of oscillations. The black hole spin barely impacts the NS oscillations in amplitude or frequency.

In [12]:
#Density Oscillations and FT Plots - Thesis figures

i=0
filename = 'MinSearch0.asc'
fig, (ax1, ax2) = plt.subplots(1,2,figsize=(15,7))
#color = [ 'k','b','orangered','k','b','orangered' ,'k','k','k','brown','m','coral','forestgreen']
#ls = ['-', '--', ':','-','-','-','--',':','-.']


for i, direc in enumerate(sorted(compare_nsbh.keys())):
    
    if compare_nsbh[direc]['type']=='BBH':
        continue
        
    labelname = compare_nsbh[direc]['label']
    datadir = os.path.join(compare_nsbh[direc]['dir'], 'Summary/Data')
    color=compare_nsbh[direc]['color']
    filepath = os.path.join(datadir, filename)
    
    t, max_rho = np.loadtxt(filepath, unpack=True, usecols=(1,5))   
    idx_stablestar = np.intersect1d(np.where(t>50), np.where(t<150))
    
    mtotal = 1.35*(1+5)
        
    idx_fft = np.intersect1d(np.where(t>0), np.where(t<150))
    
    rho_fftfreq = np.fft.fftfreq(np.size(max_rho[idx_fft]), d=(t[2]-t[1]))
    rho_fft = np.fft.fft(max_rho[idx_fft]/max_rho[0])
    
    ax1.plot(t, (max_rho)/max_rho[0], ls='-',c=color, label=labelname)
    ax1.set_xlabel('T/M')
    ax1.set_ylabel(r'$ \rho_{max}/\rho_{init}$')
    ax1.legend()
    #ax1.set_ylim(0.88, 1.5)
    ax1.set_xlim(0,150)
    ax1.grid(True)

    
     
    ax2.plot(rho_fftfreq/(mtotal/s1_in_Msun), rho_fft.real, ls='-',c=color, label=labelname)
    ax2.set_xlabel('Frequency (Hz)')
    ax2.set_ylabel(r'$FFT[\rho_{max}/\rho_{init}]$')
    ax2.set_ylim(-50,70)
    ax2.set_xlim(500,2000)
    ax2.legend()
    ax2.grid(True)

    i+=1  
plt.tight_layout() 

#plt.savefig(os.path.join(figdir, 'ReviewFigures/RevFig_DensityOscillations.png'), dpi=100)
plt.show()
plt.close()
  
    

Orbital Dynamics


Below is comparison of trajectories and orbital separation. Note that the initial momenta values differ between each case which can lead to differences in eccentricity, trajectories and number of orbits. Few key observations -

  • The merger time and number of orbits increase with spin and/or eccentricity.
  • Higher spin leads to larger eccentricity as reflected in oscillations in the orbital separation and require better eccentricity reduction techniques.
  • The BH velocity and COM offset also increases with spin of the BH at merger.
  • System facing strongest disruption ($a_{BH} = 0.75$) have highest BH speed at the merger and also suffers maximum changes in the speed. As such, its COM offset remains roughly same even after merger.
In [13]:
#Binary trajectory parameters
mpl.rc('axes', labelsize=16, grid=True)
plt.figure(figsize=(15,20))
for i, key in enumerate(sorted(compare_nsbh.keys())):
    if (compare_nsbh[key]['type']=='NSBH'):
        datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
        label_obj2='NS'
    elif(compare_nsbh[key]['type']=='BBH'):
        datadir = compare_nsbh[key]['dir']
        label_obj2='BH'
        
    t_postmerger = compare_nsbh[key]['t_merger_src']
    label = compare_nsbh[key]['label']
    bh1_search= 'ShiftTracker0.asc'
    filepath = os.path.join(datadir, bh1_search)

    if not os.path.exists(filepath): raise RuntimeError("%s file not found"%filepath)
    t1,x1,y1,z1 = np.loadtxt(filepath, unpack=True, usecols=(1,2,3,4))

    ns2_search= 'volume_integrals-GRMHD.asc'
   
    bh2_search= 'ShiftTracker1.asc'
    if (compare_nsbh[key]['type']=='NSBH'):
        filepath = os.path.join(datadir, ns2_search)
        t2,x2,y2,z2,rest_mass = np.loadtxt(filepath, unpack=True, usecols=(0,1,2,3,4))
        x2 = x2/rest_mass
        y2 = y2/rest_mass
        z2 = z2/rest_mass
       
    elif(compare_nsbh[key]['type']=='BBH'):
        filepath = os.path.join(datadir, bh2_search)
        t2,x2,y2,z2 = np.loadtxt(filepath, unpack=True, usecols=(1,2,3,4))

    plt.subplot(3,2,i+1)
    plt.plot(x1[t1<t_postmerger],y1[t1<t_postmerger], label='BH')
    plt.plot(x2[t2<t_postmerger],y2[t2<t_postmerger], label=label_obj2)
    plt.xlabel('x/M')
    plt.ylabel('y/M')
    plt.title(label, fontsize=16)
    plt.legend()
plt.tight_layout()
plt.show()
plt.close()

Final BH Offset and Coordinate velocity at merger

  • BH offset increases with spin
  • Speed increases with spin
In [14]:
#Orbital Parameters
ls=['-','--',':','-.', '--',':']
c = ['r','b','k','m','coral','forestgreen']

orbdyn_dict = {}
print('-'*110)
print("System \t \t (T_peak - R)\t \t(x, y) of BH \t \t BH Offset \t (vx, vy) of BH \t Speed")
print('-'*110)
for i, key in enumerate(sorted(compare_nsbh.keys())):
    
    if (compare_nsbh[key]['type']=='NSBH'):
        datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
        
    elif(compare_nsbh[key]['type']=='BBH'):
        datadir = compare_nsbh[key]['dir']
        
   
    label = compare_nsbh[key]['label']
    t_merger = compare_nsbh[key]['t_merger_src']
    bhsearch= 'ShiftTracker0.asc'
    filepath = os.path.join(datadir, bhsearch)

    if not os.path.exists(filepath): raise RuntimeError("%s file not found"%filepath)
    t1,x1,y1,z1,vx1,vy1,vz1 = np.loadtxt(filepath, unpack=True, usecols=(1,2,3,4,5,6,7))
        
    x1 =  x1 [t1<t_merger]
    y1 =  y1 [t1<t_merger]
    z1 =  z1 [t1<t_merger]
    vx1 = vx1[t1<t_merger]
    vy1 = vy1[t1<t_merger]
    vz1 = vz1[t1<t_merger]
       
    t1 = t1[t1<t_merger]

    t_com, x_com, y_com, z_com, vx_com, vy_com, vz_com = t1[-1], x1[-1], y1[-1], z1[-1],\
                                                        vx1[-1], vy1[-1], vz1[-1]  
    com_offset = (x_com**2 + y_com**2 + z_com**2)**0.5
    com_velocity =  np.sqrt(vx_com**2 + vy_com**2 + vz_com**2)
    print("%15s \t %.3g \t \t (%.2g, %.2g) \t%.3g \t \t(%.2g, %.2g) \t %.2g"%(label, t_com, x_com, y_com,\
                                                                com_offset,vx_com, vy_com, com_velocity))
  
    ns2_search= 'volume_integrals-GRMHD.asc'
    if (not os.path.exists(os.path.join(datadir, ns2_search))) and compare_nsbh[key]['type']=='NSBH':
        raise NameError("%s:Volume Integrals file not found. Minsearch will be used. "%compare_nsbh[key]['label'])
        
    
    bh2_search= 'ShiftTracker1.asc'
    if (compare_nsbh[key]['type']=='NSBH'):
        filepath = os.path.join(datadir, ns2_search)
        t2,x2,y2,z2,rest_mass = np.loadtxt(filepath, unpack=True, usecols=(0,1,2,3,4))
        x2 = x2[t2<t_merger]/rest_mass[t2<t_merger]
        y2 = y2[t2<t_merger]/rest_mass[t2<t_merger]
        z2 = z2[t2<t_merger]/rest_mass[t2<t_merger]
        t2 = t2[t2<t_merger]
        vx2 = np.gradient(x2)/np.gradient(t2)
        vy2 = np.gradient(y2)/np.gradient(t2)
        vz2 = np.gradient(z2)/np.gradient(t2)
            
       
    elif (compare_nsbh[key]['type']=='BBH'):
        filepath = os.path.join(datadir, bh2_search)
        t2,x2,y2,z2,vx2,vy2,vz2 = np.loadtxt(filepath, unpack=True, usecols=(1,2,3,4,5,6,7))
        x2 =  x2 [t2<t_merger]
        y2 =  y2 [t2<t_merger]
        z2 =  z2 [t2<t_merger]
        vx2 = vx2[t2<t_merger]
        vy2 = vy2[t2<t_merger]
        vz2 = vz2[t2<t_merger]
        
        t2 = t2[t2<t_merger]
        
    t_common = np.intersect1d(t1, t2)

    idx_st = np.isin(t1,t_common)
    idx_ms = np.isin(t2,t_common)

    sep = np.sqrt((x1[idx_st]-x2[idx_ms])**2 + (-y2[idx_ms]+y1[idx_st])**2 + (-z2[idx_ms]+z1[idx_st])**2)
    orbphi = np.unwrap(np.arctan2(-y2[idx_ms]+y1[idx_st],-x2[idx_ms]+x1[idx_st]))
    orbtheta = np.arccos(np.divide(z1[idx_st]-z2[idx_ms],sep))

    vx = vx1[idx_st]-vx2[idx_ms]
    vy = vy1[idx_st]-vy2[idx_ms]
    vz = vz1[idx_st]-vz2[idx_ms]
    
    #Velocities computed for BH are much more accurate as it does not use FD and computed from Shift. As we lack 
    #data for every iteration in NS + use second order FD, velocities tend to have higher numerical errors. Hence
    #compute rdot and phidot directly from numerical derivative of the functions
    
    #rdot = (vx*np.cos(orbphi) + vy*np.sin(orbphi))*np.sin(orbtheta) + vz*np.cos(orbtheta)
    #phidot =  np.divide((vy*np.cos(orbphi) - vx*np.sin(orbphi)), (sep*np.sin(orbtheta)))  
    
    rdot = np.gradient(sep)/np.gradient(t1[idx_st])
    phidot = np.gradient(orbphi)/np.gradient(t1[idx_st])
    
    orbdyn_dict[key] = {}
    orbdyn_dict[key]['time'] = t_common
    orbdyn_dict[key]['separation'] = sep
    orbdyn_dict[key]['phase'] = orbphi
    orbdyn_dict[key]['rdot'] = rdot
    orbdyn_dict[key]['omega'] = phidot
   
--------------------------------------------------------------------------------------------------------------
System 	 	 (T_peak - R)	 	(x, y) of BH 	 	 BH Offset 	 (vx, vy) of BH 	 Speed
--------------------------------------------------------------------------------------------------------------
$a_{BH} = -0.5$ 	 156 	 	 (-0.038, 0.0013) 	0.0378 	 	(-0.0032, 0.0013) 	 0.0035
   $a_{BH} = 0$ 	 329 	 	 (-0.12, -0.023) 	0.127 	 	(0.0024, 0.0027) 	 0.0036
$a_{BH} = 0.25$ 	 453 	 	 (-0.17, -0.066) 	0.182 	 	(0.0021, 0.0035) 	 0.0041
 $a_{BH} = 0.5$ 	 584 	 	 (-0.22, -0.21) 	0.306 	 	(0.006, 0.011) 	 0.013
$a_{BH} = 0.75$ 	 714 	 	 (-0.12, -0.67) 	0.68 	 	(0.014, 0.029) 	 0.032

Final BH Offset and Coordinate velocity at 300M after merger

  • BH Offset continues to increase but the differences with offset at merger decrease with spin
  • Coordinate speed drastically decreases for large spin values and do not follow the same trend as before.
In [15]:
#Orbital Parameters
ls=['-','--',':','-.', '--',':']
c = ['r','b','k','m','coral','forestgreen']

sep_list = []
orbphase_list = []
rdot_list = []
phidot_list = []
t_common_list = []
label_list = []

print('-'*110)
print("System \t \t (T_peak - R)+300 \t (x, y) of BH \t \tBH Offset \t (vx, vy) of BH \t Speed")
print("-"*110)
for i, key in enumerate(sorted(compare_nsbh.keys())):
    
    if (compare_nsbh[key]['type']=='NSBH'):
        datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
        t_postmerger = compare_nsbh[key]['t_merger_src']
    elif(compare_nsbh[key]['type']=='BBH'):
        datadir = compare_nsbh[key]['dir']
        t_postmerger = compare_nsbh[key]['t_merger_src']
   
    label = compare_nsbh[key]['label']
    bhsearch= 'ShiftTracker0.asc'
    filepath = os.path.join(datadir, bhsearch)

    if not os.path.exists(filepath): raise RuntimeError("%s file not found"%filepath)
    t1,x1,y1,z1,vx1,vy1,vz1 = np.loadtxt(filepath, unpack=True, usecols=(1,2,3,4,5,6,7))
    
    if t1[-1]>t_postmerger:
        idx_300 = np.amin(np.where(t1>=t_postmerger+300))
        t_com, x_com,y_com, z_com, vx_com, vy_com, vz_com = t1[idx_300], x1[idx_300], y1[idx_300], z1[idx_300],\
                                                            vx1[idx_300], vy1[idx_300], vz1[idx_300]  
        com_offset = (x_com**2 + y_com**2 + z_com**2)**0.5
        com_velocity = np.sqrt(vx_com**2 + vy_com**2 + vz_com**2)
        if y_com<0:
            print("%15s \t %.4g \t \t  (%.2f, %.2f) \t%.3g \t\t (%.2g, %.2g) \t %.2g"%(compare_nsbh[key]['label'],\
                                                                                    t_com, x_com, y_com, com_offset, \
                                                                                    vx_com, vy_com, com_velocity ))
        else:
            print("%15s \t %.4g \t \t  (%.2f, %.2f) \t%.3g \t\t (%.2g, %.2g) \t %.2g"%(compare_nsbh[key]['label'],\
                                                                                    t_com, x_com, y_com,  com_offset, \
                                                                                        vx_com, vy_com, com_velocity ))
       
   
--------------------------------------------------------------------------------------------------------------
System 	 	 (T_peak - R)+300 	 (x, y) of BH 	 	BH Offset 	 (vx, vy) of BH 	 Speed
--------------------------------------------------------------------------------------------------------------
$a_{BH} = -0.5$ 	 455.6 	 	  (-0.10, 0.12) 	0.152 		 (-0.00013, 0.00042) 	 0.00044
   $a_{BH} = 0$ 	 629.3 	 	  (-0.33, -0.08) 	0.343 		 (-0.00074, -0.00024) 	 0.00078
$a_{BH} = 0.25$ 	 753.3 	 	  (-0.31, -0.17) 	0.353 		 (-0.00053, -0.00042) 	 0.00068
 $a_{BH} = 0.5$ 	 883.9 	 	  (-0.30, -0.27) 	0.402 		 (-0.00029, -0.00049) 	 0.00057
$a_{BH} = 0.75$ 	 1014 	 	  (-0.26, -0.64) 	0.693 		 (7.8e-05, -0.00071) 	 0.00071

Orbital Separation and Phase

In [16]:
#Orbital Separation and Phase - NSBH vs BBH
fig, ax = plt.subplots(2,1,figsize=(10,10))

print('-'*100)
print("Simulation \t Merger Time \t Number of Orbits")
print('-'*100)
keys = compare_nsbh.keys()
for i,key in enumerate(sorted(compare_nsbh.keys())):
    
    label = compare_nsbh[key]['label']
    color = compare_nsbh[key]['color']
    ls    = compare_nsbh[key]['ls']
        
        
    print("{:<16s} \t {:>4.3g} \t \t{:>4.3g} ".format(label,orbdyn_dict[key]['time'][-1], orbdyn_dict[key]['phase'][-1]/2./np.pi))
    
    #Derivatives of mixed binary can be extremely noise. So instead of using the direct data, we perform 
    #fitting using spline interpolation and use them to compute the derivatives
    sep_splinefit = splrep(orbdyn_dict[key]['time'],orbdyn_dict[key]['separation'],k=5,s=3)
    phase_splinefit = splrep(orbdyn_dict[key]['time'], orbdyn_dict[key]['phase'],k=5,s=3)
    
    orbdyn_dict[key]['sepfit'] = sep_splinefit
    orbdyn_dict[key]['phasefit'] = phase_splinefit
    
    
    if 'nsbh' in key: 
        #td_radius = compare_nsbh[key]['tidal_radius']
        #Compute time of reaching tidal radius
        time_key = orbdyn_dict[key]['time']
        sep_key = orbdyn_dict[key]['separation']
        #time_tidalrad = time_key[np.amin(np.argwhere(sep_key<=td_radius))]
   
    ax[0].plot(orbdyn_dict[key]['time'], orbdyn_dict[key]['separation'], ls='-',c=color,label=label)
    #plt.plot(orbdyn_dict[key]['time'], splev(orbdyn_dict[key]['time'], sep_splinefit), ls='--',c='k', lw=1,label='Spline Fit')
    ax[0].set_xlabel('T/M', fontsize=16)
    ax[0].set_ylabel('Separation', fontsize=16)
    ax[0].legend()
    
    ax[1].plot(orbdyn_dict[key]['time'], orbdyn_dict[key]['phase'], ls='-',c=color,label=label)
    #plt.plot(orbdyn_dict[key]['time'], splev(orbdyn_dict[key]['time'], phase_splinefit), ls=':',c='k',label='Spline Fit')
    ax[1].set_xlabel('T/M', fontsize=16)
    ax[1].set_ylabel('Orbital Phase', fontsize=16)
    ax[1].legend()
   
print('\n\n')
plt.tight_layout()
fig.suptitle('Figure 2.1 - NSBH vs BBH', fontsize=20, y=1.02)
plt.show()
plt.close()
----------------------------------------------------------------------------------------------------
Simulation 	 Merger Time 	 Number of Orbits
----------------------------------------------------------------------------------------------------
$a_{BH} = -0.5$  	  155 	 	1.79 
$a_{BH} = 0$     	  329 	 	3.56 
$a_{BH} = 0.25$  	  453 	 	4.72 
$a_{BH} = 0.5$   	  584 	 	 5.9 
$a_{BH} = 0.75$  	  714 	 	6.95 



Psi4 - (2,2) Mode Comparison

In [17]:
#Re(Psi4) plots: 
import matplotlib.ticker as plticker

fig, ax = plt.subplots(1, 1, figsize=(15,7))
 

for i, key in enumerate(sorted(compare_nsbh.keys())):
    
    labelname = compare_nsbh[key]['label']
    color=compare_nsbh[key]['color']
    ls = compare_nsbh[key]['ls']
    
    if (compare_nsbh[key]['type']=='NSBH'):
        datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
    elif(compare_nsbh[key]['type']=='BBH'):
        datadir = compare_nsbh[key]['dir']

    filepath = os.path.join(datadir, 'Ylm_WEYLSCAL4::Psi4r_l2_m2_r70.00.asc')
    if not os.path.exists(filepath):
        filepath = os.path.join(datadir, 'mp_Psi4r_l2_m2_r70.00.asc')
   
    t22, psi4_re22, psi4_im22 = np.loadtxt(filepath, unpack=True, usecols=(0,1,2))
    amp = np.sqrt(psi4_re22**2 + psi4_im22**2)
    r = float(((filepath.split('.asc')[0]).split('_')[-1]).split('r')[-1])
    
    
    ax.plot(t22, r*psi4_re22,c=color,ls='-', label=labelname)
    #ax.plot(t22,amp,ls=ls,c=color)
    ax.set_xlabel(r'$T/M$')
    ax.set_ylabel(r'$Re(\psi_4)\,r\,M$')
    ax.set_xlim(0,900)
    ax.legend(loc=3)
    ax.grid(True)
    ax.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
    loc = plticker.MultipleLocator(base=100) # this locator puts ticks at regular intervals
    ax.xaxis.set_major_locator(loc)
     

plt.tight_layout()
#plt.savefig(os.path.join(figdir, 'C3_Psi4Plots.png'), dpi=200)
plt.show()
plt.close()

Psi4 Amplitude

In [18]:
#Psi4 Amplitude plots: BBH vs NSBH

i=0


color = ['k','b', 'orangered','r','m','coral','forestgreen']
ls = ['-', '--', ':','-.','--',':']

fig=plt.figure(figsize=(15,8))
for i, key in enumerate(sorted(compare_nsbh.keys())):
    
    labelname = compare_nsbh[key]['label']
    col = compare_nsbh[key]['color']
    
    if (compare_nsbh[key]['type']=='NSBH'):
        datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
        ls = '-'
    elif(compare_nsbh[key]['type']=='BBH'):
        datadir = compare_nsbh[key]['dir']
        ls=':'
    
    filepath = os.path.join(datadir, 'Ylm_WEYLSCAL4::Psi4r_l2_m2_r70.00.asc')
    if not os.path.exists(filepath):
        filepath = os.path.join(datadir, 'mp_Psi4r_l2_m2_r70.00.asc')
    
    R = float(((filepath.split('.asc')[0]).split('_')[-1]).split('r')[1])
    t22, psi4_re22, psi4_im22 = np.loadtxt(filepath, unpack=True, usecols=(0,1,2))
    psi4_re22 = psi4_re22[t22>100]
    psi4_im22 = psi4_im22[t22>100]
    t22       = t22[t22>100]   
    
    amp = np.sqrt(psi4_re22**2 + psi4_im22**2)
    max_amp_idx = np.argmax(amp==np.amax(amp))
    
    plt.subplot(121)
    plt.plot(t22-t22[max_amp_idx], R*amp,ls=ls, c=col, label=labelname)
    
    plt.xlabel(r'$T - R_{ex}(M)$')
    plt.ylabel(r'$r \, |\Psi_4|M$')
    plt.xlim(-50,100)
    plt.ylim(1e-5, 0.03)
    plt.legend()
    plt.grid(True)
    
    if "-0.5" in labelname:
        continue
    plt.subplot(122)
    plt.plot(t22, R*amp,ls=ls, c=col, label=labelname)
    
    plt.xlabel(r'$T/M$')
    plt.ylabel(r'$r \, |\Psi_4|M$')
    plt.xlim(150,300)
    plt.ylim(0,0.002)
    plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
   
    plt.grid(True)
    plt.legend()
    
plt.tight_layout()

#fig.suptitle('Figure 2: BBH vs NSBH Psi4 (2,2) Mode', fontsize=16, y=1.08)
#plt.savefig(os.path.join(figdir,'C3_Psi4Amp.png'), dpi=100)
plt.show()
plt.close()
Adding an axes using the same arguments as a previous axes currently reuses the earlier instance.  In a future version, a new instance will always be created and returned.  Meanwhile, this warning can be suppressed, and the future behavior ensured, by passing a unique label to each axes instance.

Quasinormal Fits

Each column represents the $\psi_4$ mode as a function of time and rows look at log(amplitude) and phase. The dashed lines correspond to the simulation data where colors indicate the simulation. Solid lines are fit to the data (not quasinormal ringdown fits!).

To Check - If (3,2) mode in BBH is also equally noisy.

In [19]:
#Computing the quasinormal frequencies with given mass and spins using Ringdown fits
import pandas as pd
#Reference - PHYSICAL REVIEW D73,064030 (200

def fit_omega(a_bh, M_bh, f1, f2, f3):
    
    omega_lmn = (f1 + f2*(1.-a_bh)**f3)/M_bh
    return omega_lmn

def fit_tau(a_bh, omega_lmn, q1, q2, q3):
    tau_lmn = (2./omega_lmn)*(q1+q2*(1.-a_bh)**q3)
    return tau_lmn

    
f1_dict = {'220':1.5251, '210':0.6,'200':0.4437, '2-10':0.3441, '2-20':0.2938,\
           '330':1.8956, '320':1.1481,'310':0.8345,'300':0.6873,'3-10':0.5751,'3-20':0.5158,'3-30':0.4673}
f2_dict = {'220':-1.1568, '210':-0.2339,'200':-0.0739, '2-10':0.0293, '2-20':0.0782,\
           '330':-1.3043, '320':-0.5552,'310':-0.24205,'300':-0.09282,'3-10':0.02508,'3-20':0.8195,'3-30':0.1296}
f3_dict = {'220':0.1292, '210':0.4175,'200':0.335, '2-10':2.001, '2-20':1.3546,\
           '330':0.1818, '320':0.3002,'310':0.4095,'300':0.3479,'3-10':3.136,'3-20':1.4084,'3-30':1.3255}

q1_dict = {'220':0.7, '210':-0.3,'200':4.0, '2-10':2.0, '2-20':1.67,
           '330':0.9, '320':0.8313,'310':23.845,'300':6.7841,'3-10':3.0464,'3-20':2.9,'3-30':2.55}
q2_dict = {'220':1.4187, '210':2.3561,'200':-1.955, '2-10':0.1078, '2-20':0.4192,
           '330':2.343, '320':2.3773,'310':-20.724,'300':-3.6112,'3-10':0.1162,'3-20':0.3356,'3-30':0.6576}
q3_dict = {'220':-0.499, '210':-0.2277,'200':0.142, '2-10':5.0069, '2-20':1.47,
           '330':-0.481, '320':-0.3655,'310':0.03837,'300':0.0948,'3-10':-0.2812,'3-20':2.305,'3-30':1.3378}

#From results obtained below
q2_bbh  = ['BBH', 2.0,0.6137, 0.9591]
q2_nsbh = ['NSBH',2.0,0.6846, 0.9672]
q3_bbh  = ['BBH', 3.0,0.5405, 0.9712]
q3_nsbh = ['NSBH',3.0,0.564, 0.9746]
q5_bbh  = ['BBH', 5.0,0.4166, 0.9824]
q5_nsbh = ['NSBH',5.0,0.4204, 0.9834]


remnant_data = pd.DataFrame(np.array((q2_bbh,  q2_nsbh, q3_bbh, q3_nsbh, q5_bbh, q5_nsbh)),columns=('System','q','a','Mh') )


for key in ['210', '220','330']:#,'440']:
    f1 = f1_dict[key]
    f2 = f2_dict[key]
    f3 = f3_dict[key]
    
    q1 = q1_dict[key]
    q2 = q2_dict[key]
    q3 = q3_dict[key]
    
    a_bh = pd.to_numeric(remnant_data['a'])
    M_bh = pd.to_numeric(remnant_data['Mh'])
    
    
    omega = fit_omega(a_bh, M_bh, f1, f2, f3)
    tau = fit_tau(a_bh, omega, q1, q2, q3)
    
    remnant_data['Omega_%s'%key] = omega
    remnant_data['tau_%s'%key] = tau
    #remnant_data['MOmega_%s'%key] = M_bh*omega
    #remnant_data['Tau_by_M%s'%key] = M_bh*1./tau
    
    
    
In [20]:
#QNM Fitting - Thesis Figure 
#Psi4 plots: BBH vs NSBH
from collections import defaultdict

i=0


color = ['k','b', 'orangered','m','forestgreen','r', 'b', 'k','coral']
ls = ['-', '--', ':','-.','--',':']


lm_modes = {'2,2':(2,2), '2,1':(2,1),'3,3':(3,3),'3,2':(3,2), '4,4':(4,4)}
y_amp_max = [1e-4, 5e-4, 1e-4, 1e-4]
y_phase_max = [150, 200,250,250]
y_phase_min = [-25, -50,-20,50]

qnm_tau_dict = defaultdict(dict)
qnm_freq_dict = defaultdict(dict)

fig = plt.figure(figsize=(14,20))

nsbh_iter = 0
  
for lm_idx, lm_key in enumerate(sorted(lm_modes.keys())):

    lm=lm_modes[lm_key]
   
    for i, key in enumerate(sorted(compare_nsbh.keys())):
   
        labelname = compare_nsbh[key]['label']
        col=compare_nsbh[key]['color']#color[i]
        
        if (compare_nsbh[key]['type']=='NSBH'):
            datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
            nsbh_iter+=1
            #remnant_data_nsbh = remnant_data[remnant_data['System']=='NSBH']
 
          
        q = 5
        #q_arr = pd.to_numeric(remnant_data_nsbh['q'])
        #remnant_data_nsbh = remnant_data_nsbh[q_arr==q]
        
        #omega = np.float(remnant_data_nsbh['Omega_%d%d0'%(lm[0],lm[1])])
        #tau =  np.float(remnant_data_nsbh['tau_%d%d0'%(lm[0],lm[1])])
       
        filepath = os.path.join(datadir, 'Ylm_WEYLSCAL4::Psi4r_l%d_m%d_r70.00.asc'%(lm[0],lm[1]))
        if not os.path.exists(filepath):
            filepath = os.path.join(datadir, 'mp_Psi4r_l%d_m%d_r70.00.asc'%(lm[0],lm[1]))
        
        t22, psi4_re22, psi4_im22 = np.loadtxt(filepath, unpack=True, usecols=(0,1,2))
    
        amp = np.sqrt(psi4_re22**2 + psi4_im22**2)
        phase = -1*np.unwrap(np.arctan2(psi4_im22,psi4_re22))
        
        #Compute QNM Modes
        #For QNM - A = constant*exp(-t\Ï„) and phase = omega*t+phi_0
        
        #Find the fitting time interval
        max_amp_idx = np.argmax(amp==np.amax(amp))
        qnm_begin_idx = np.amin(np.where(t22 > t22[max_amp_idx]+30))
        qnm_end_idx = np.amin(np.where(t22 > t22[max_amp_idx]+50))
        
        
        #Define the fitting region
        qnm_time_interval     = t22[qnm_begin_idx:qnm_end_idx]
        qnm_amp_interval      = amp[qnm_begin_idx:qnm_end_idx]
        qnm_phase_interval    = phase[qnm_begin_idx:qnm_end_idx]
         
        #Create the fit for amplitude and phase
        logamp_fit_coeff   = np.polyfit(qnm_time_interval, np.log(qnm_amp_interval),1)
        logamp_fit_data   = np.poly1d(logamp_fit_coeff)
        #logamp_fit_berti  = lambda x: -(1./tau)*x + (np.log(amp)[qnm_begin_idx] + (1./tau)*t22[qnm_begin_idx])
        
        
        phase_fit_coeff    = np.polyfit(qnm_time_interval, qnm_phase_interval,1)
        phase_fit_data  = np.poly1d(np.array([phase_fit_coeff[0], phase_fit_coeff[1]]))#np.poly1d(phase_fit_coeff)
        #phase_fit_berti  = lambda x: omega*x + (phase[qnm_begin_idx] - omega*t22[qnm_begin_idx])
        
        #Compute the Quasinormal Frequency and Decay Time constant from the fits
        logamp_fit_slope = logamp_fit_coeff[0]
        
        qnm_tau = -1./logamp_fit_slope
        qnm_freq = phase_fit_coeff[0]/2./np.pi
        
        qnm_tau_dict[lm_key][labelname] = qnm_tau
        qnm_freq_dict[lm_key][labelname] = qnm_freq
        
        #Plot the data and the fits
        plt.subplot(5,2,2*lm_idx+1)
        amp_plot = plt.semilogy(t22-t22[max_amp_idx], amp, ls='--', c=col, label=labelname)
        
        plt.semilogy(qnm_time_interval-t22[max_amp_idx], np.exp(logamp_fit_data(qnm_time_interval)),ls='-', c=col)
        #plt.semilogy(qnm_t_int-t22[max_amp_idx], np.exp(logamp_fit_berti(qnm_t_int)),ls='-', c=color[lm_idx])
    
        if nsbh_iter==3:plt.xlabel(r'$T/M$')
        plt.ylabel(r'$|\psi^{ %d,%d}|$'%(lm[0],lm[1]))
        plt.xlim(0,130)
        plt.ylim(1e-9,4e-4)
        plt.legend()
        plt.grid(True)
        
        
        plt.subplot(5,2,2*lm_idx+2)
        plt.plot(t22 - t22[max_amp_idx], phase,ls='--', c=col, label=labelname+'(%d,%d)'%(lm[0],lm[1]))
        plt.plot(qnm_time_interval - t22[max_amp_idx], phase_fit_data(qnm_time_interval), ls='-', c=col)#, label=labelname+'Fit')
        
        #plt.plot(qnm_t_int - t22[max_amp_idx], phase_fit_berti(qnm_t_int), ls='-', c=color[lm_idx])#, label=labelname+'Fit')
        
        if lm_idx==2:plt.xlabel(r'$T/M$')
        plt.ylabel(r'$\phi^{ %d,%d}$'%(lm[0],lm[1]))
        plt.xlim(0,130)
        #plt.ylim(y_phase_min[nsbh_iter-1], y_phase_max[nsbh_iter-1])
        plt.legend(loc=2)
        plt.grid(True)
    
#fig.suptitle('(%s) Mode'%lm_key, fontsize=16, y=1.08)
plt.tight_layout()
#plt.savefig(os.path.join(figdir, 'QNMFitting_NSBH.png'), dpi=100)
plt.show()
plt.close()
In [ ]:
# QNM Tau & QNM Frequency Comparison - Thesis Figure

color = ['C0','C1','C2','C3','C4','C5']
fig=plt.figure(figsize=(15,6))
plt.subplot(121)
for key1 in sorted(qnm_tau_dict.keys()):
    print("\n\nMode: l=%s, m=%s \n"%(key1[0], key1[2]))
    print("Simulation \t QNM Tau \t QNM Freqeuncy")
    print("-"*50)

    for i,key2 in enumerate(sorted(qnm_tau_dict[key1].keys())):

        print("%10s \t %g \t %g "%(key2, qnm_tau_dict[key1][key2], qnm_freq_dict[key1][key2]*2*np.pi))
        if key1=='2,2':label=key2
        else: label= "_nolegend_"
                
        if key2=='Q2-NSBH':continue
        plot = plt.scatter(key1, qnm_tau_dict[key1][key2], c=color[i],label=label)
    
plt.xlabel('(l,m) Mode')
plt.ylabel(r'$\tau_{QNM}$')
plt.grid(True)
plt.legend()

color = ['C0','C1','C2','C3','C4','C5']
plt.subplot(122)
for key1 in sorted(qnm_freq_dict.keys()):
    for i,key2 in enumerate(sorted(qnm_freq_dict[key1].keys())):
        
        if key2=='Q2-NSBH':continue
        if key1=='2,2':label=key2
        else: label= "_nolegend_"
        plot = plt.scatter(key1, qnm_freq_dict[key1][key2], c=color[i],label=label)
    
print('\n\n')
plt.xlabel('(l,m) Mode')
plt.ylabel(r'$f_{QNM}$')
plt.grid(True)
plt.legend()
plt.tight_layout()
#fig.suptitle('Figure 7: QNM Quatities', y=1.08, fontsize=16)
plt.savefig(os.path.join(figdir,'QNMModes.png'), dpi=100)
#plt.show()
plt.close()

Radiated Energy, Angular Momentum and BH Kicks

Radiated Quantities


Three key observations -

  1. Radiated energy and angular momentum after merger increase with BH spin.
  2. At the merger, BH kicks are highest for non-spinning case and gradually decrease for each orientation of the spin. In contrast, post merger kicks are highest for anti-aligned spin.
  3. BH suffers a strong retardation for very high spin in +z direction which is reflected in decrease in the speed of the black hole. In contrast the speeds increase for the anti-aligned case!
In [21]:
#Radiated Quantities

from scipy.constants import c, G
print('-'*100)
print('Simulation \t  Time after merger \t  Radiated Energy \t Luminosity \t \tRadiated Jz  ')
print('-'*100)
for i, key in enumerate(sorted(compare_nsbh.keys())): #[0:3:2])):
    
    labelname = compare_nsbh[key]['label']

    if (compare_nsbh[key]['type']=='NSBH'):
        t_merge = compare_nsbh[key]['t_merger_src']
        datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
    elif(compare_nsbh[key]['type']=='BBH'):
        datadir = compare_nsbh[key]['dir']
        t_merge = compare_nsbh[key]['t_merger_src']
    
    filepath = os.path.join(datadir, 'psi4radAway_l6_r70.00.asc')
    
    #filepath = os.path.join(datadir, 'psi4analysis_r70.00.asc')
    if not os.path.exists(filepath):continue
    t, Erad, Kx, Ky, Kz, Edot, Pxdot, Pydot, Pzdot, Jxdot, Jydot, Jzdot, Jx, Jy, Jz\
             = np.loadtxt(filepath, unpack=True, usecols=(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14))
    time_centered = t - (t_merge + 70)
    try:
        idx = np.amin(np.where(time_centered>200))
    except ValueError as err:
        idx=-1
    
    
    print ('{:<16s} \t {:>4.3g} \t \t {:>4.3g} \t \t{:>4.3g} \t {:>12g} '.format(labelname, time_centered[idx], Erad[idx], Edot[idx], Jz[idx]))
    
   
print('\n\n')
print('-'*100)
print('Simulation \t Merger Time  \tBH Mass \t Kick (km/s) \t Merger Time+200 BH Mass \t Kick  (km/s)')
print('-'*100)
for i, key in enumerate(sorted(compare_nsbh.keys())): #[0:3:2])):
    
    labelname = compare_nsbh[key]['label']
    if (compare_nsbh[key]['type']=='NSBH'):
        datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
        rad_arr = np.array((50,60,70,80,90,100,110,120,130,140))
        t_merge = compare_nsbh[key]['t_merger_src']
       
    elif(compare_nsbh[key]['type']=='BBH'):
        datadir = compare_nsbh[key]['dir']
        rad_arr = np.array((30,40,50,60,70,75,80,90,100,115,130,145,160,175,190,200))
        t_merge = compare_nsbh[key]['t_merger_src']
       
    
    #filepath = os.path.join(datadir, 'psi4analysis_r70.00.asc')
    filepath = os.path.join(datadir, 'psi4radAway_l6_r70.00.asc')
    if not os.path.exists(filepath):continue
    t, Erad, Kx, Ky, Kz, Edot, Pxdot, Pydot, Pzdot, Jxdot, Jydot, Jzdot, Jx, Jy, Jz\
             = np.loadtxt(filepath, unpack=True, usecols=(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14))
    Kmag = (Kx**2 + Ky**2 + Kz**2)**0.5
    time_centered = t - (t_merge+70)
    
    idx_merger = np.amin(np.where(time_centered>0))
    idx_postmerger = np.amin(np.where(time_centered>200))
    
    bh_diag_file = os.path.join(datadir, 'BH_diagnostics.ah1.gp')
    ihspin_file = os.path.join(datadir, 'ihspin_hn_0.asc')
    t_ah, irr_mass = np.loadtxt(bh_diag_file, unpack=True, usecols=(1,26))
    t_ih, sx, sy, sz = np.loadtxt(ihspin_file, unpack=True, usecols=(0,1,2,3))
    
    
    ah_idx_merger = np.amin(np.where(t_ah>t_merge + 0))
    ah_idx_postmerger = np.amin(np.where(t_ah>t_merge + 200))
    ih_idx_merger = np.amin(np.where(t_ih>t_merge + 0))
    ih_idx_postmerger = np.amin(np.where(t_ih>t_merge + 200))
    
    spin_mag_merger = np.linalg.norm(np.array((sx, sy, sz)).T, axis=1)[ih_idx_merger]
    hzn_mass_merger = np.sqrt(irr_mass[ah_idx_merger]**2 + 0.25*spin_mag_merger**2/irr_mass[ah_idx_merger]**2)
    
    spin_mag_postmerger = np.linalg.norm(np.array((sx, sy, sz)).T, axis=1)[ih_idx_postmerger]
    hzn_mass_postmerger = np.sqrt(irr_mass[ah_idx_postmerger]**2 + 0.25*spin_mag_postmerger**2/irr_mass[ah_idx_postmerger]**2)
   
    #print ('{:<16s}   {:>1.2g} \t {:>1.2g} \t {:>1.2g} \t {:>1.2g} \t {:>1.4g} \t {:>1.4g}'.format(labelname,t[idx],Kx[idx]/hzn_mass, Ky[idx]/hzn_mass, \
    #                                                                                   Kz[idx]/hzn_mass, hzn_mass, (Kmag[idx]/hzn_mass)*c/1000 ))
    
    print ('{:<16s}   {:>1.1f} \t \t {:>1.4g} \t {:>1.3g}  \t\t {:>1.1f}  \t {:>1.4g}  \t {:>1.4g} '.format(labelname,time_centered[idx_merger], \
                                                                hzn_mass_merger, \
                                                                (Kmag[idx_merger]/hzn_mass_merger)*c/1000,\
                                                                time_centered[idx_postmerger],hzn_mass_postmerger, \
                                                                (Kmag[idx_postmerger]/hzn_mass_postmerger)*c/1000 ))
    
----------------------------------------------------------------------------------------------------
Simulation 	  Time after merger 	  Radiated Energy 	 Luminosity 	 	Radiated Jz  
----------------------------------------------------------------------------------------------------
$a_{BH} = -0.5$  	  200 	 	 0.00676 	 	1.87e-10 	   -0.0550821 
$a_{BH} = 0$     	  200 	 	 0.00947 	 	9.37e-11 	   -0.0845095 
$a_{BH} = 0.25$  	  200 	 	 0.0111 	 	1.44e-10 	    -0.102625 
$a_{BH} = 0.5$   	  200 	 	 0.012 	 	5.69e-11 	    -0.118622 
$a_{BH} = 0.75$  	  200 	 	 0.0125 	 	4.76e-11 	     -0.13151 



----------------------------------------------------------------------------------------------------
Simulation 	 Merger Time  	BH Mass 	 Kick (km/s) 	 Merger Time+200 BH Mass 	 Kick  (km/s)
----------------------------------------------------------------------------------------------------
$a_{BH} = -0.5$    0.3 	 	 0.9868 	 87.4  		 200.1  	 0.9875  	 164.4 
$a_{BH} = 0$       0.3 	 	 0.9818 	 115  		 200.1  	 0.9841  	 70.04 
$a_{BH} = 0.25$    0.3 	 	 0.9716 	 113  		 200.1  	 0.9803  	 25.85 
$a_{BH} = 0.5$     0.3 	 	 0.9383 	 95.2  		 200.1  	 0.9719  	 35.55 
$a_{BH} = 0.75$    0.3 	 	 0.8677 	 83.5  		 200.1  	 0.9531  	 37.1 

Radiated Energy, Luminosity and Angular Momentum

  • System with highest initial BH spin emits maximum energy in gravitational radiation and follows a decreasing trend with spin. In contrast, the luminosity peaks at $a_{BH}=0.25$ and decreases for all other values. This trend differs slightly from $\psi_4^{22}$.

  • Angular momentum shows two interesting trends - first for non-negative spinning systems, there exist additional angular momentum trapped in the remaining matter. Second, the gravitational radiation is emitted for longer duration in situations with minimal stellar disruptions. In contrast, the loss of angular momentum in gravitational waves stops close to the merger for high spinning cases.

In [22]:
#Energy and Angular Momentum Plot:

i=0

plt.figure(figsize=(15,12))

color = [ 'k','b','orangered','m','forestgreen']
ls = ['-', '--', ':','-.','--',':']


for i, key in enumerate(sorted(compare_nsbh.keys())): #[0:3:2])):
    
    
    labelname = compare_nsbh[key]['label']
    col = compare_nsbh[key]['color']
    if (compare_nsbh[key]['type']=='NSBH'):
        datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
        ls='-'
        
    elif(compare_nsbh[key]['type']=='BBH'):
        datadir = compare_nsbh[key]['dir']
        ls='--'
        
    filepath = os.path.join(datadir, 'psi4radAway_l6_r70.00.asc')#'psi4radAway_l6_r70.00.asc')
    
    if not os.path.exists(filepath):
        print('File note found for {}'.format(labelname))
        
    t, Erad, Kx, Ky, Kz, Edot, Pxdot, Pydot, Pzdot, Jxdot, Jydot, Jzdot, Jx, Jy, Jz\
             = np.loadtxt(filepath, unpack=True, usecols=(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14))

    
    psi4_filepath = os.path.join(datadir, 'Ylm_WEYLSCAL4::Psi4r_l2_m2_r70.00.asc')
    if not os.path.exists(psi4_filepath):
        psi4_filepath = os.path.join(datadir, 'mp_Psi4r_l2_m2_r70.00.asc')
        
    t22, psi4_re22, psi4_im22 = np.loadtxt(psi4_filepath, unpack=True, usecols=(0,1,2))
    amp = np.sqrt(psi4_re22**2 + psi4_im22**2)
    max_amp_idx = np.argmax(amp==np.amax(amp))
    
    plt.subplot(221)
    plt.plot(t-t22[max_amp_idx], Erad,ls=ls, c=col, label=labelname)
    plt.xlabel(r'$T/M$')
    plt.ylabel(r'$E_{rad}/M$')
    plt.xlim(-400,200)
    plt.legend()
    
    plt.subplot(222)
    plt.plot(t, Edot,ls=ls, c=col, label=labelname)
    plt.xlabel(r'$T/M$')
    plt.ylabel(r'$\dot{E}_{rad}$')
    #plt.xlim(-150,100)
    plt.legend()

    plt.subplot(223)
    plt.plot(t-t22[max_amp_idx], Jz,ls=ls, c=col, label=labelname)
    plt.xlabel(r'$T/M$')
    plt.ylabel(r'$J_{rad}/M$')
    plt.xlim(-400,200)
    plt.legend()
    
    plt.subplot(224)
    plt.plot(t, Jzdot,ls=ls, c=col, label=labelname)
    plt.xlabel(r'$T/M$')
    plt.ylabel(r'$\dot{J}_{rad}$')
    #plt.xlim(-150,100)
    plt.legend()

    
plt.tight_layout()

#plt.savefig(os.path.join(figdir,'ERad.png'), dpi=100)
plt.show()
plt.close()
In [23]:
#Angular Momenta and BH Spin
import matplotlib.ticker as plticker


i=0

fig, ax = plt.subplots(2,3, figsize=(16,10))
ax = np.concatenate(ax)
color = ['r', 'b', 'k','m','coral','forestgreen']
ls = ['--', '-', ':','-.','--',':']


for i, key in enumerate(sorted(compare_nsbh.keys())):
    labelname = compare_nsbh[key]['label']
    mass_ratio = 5
    spin_bh = float(labelname.split('=')[1].split('$')[0])
    ls = compare_nsbh[key]['ls']
    color=compare_nsbh[key]['color']
    
    if spin_bh==-0.5: 
        ax1 = ax[0]
        tmax = 300
        #ymin,ymax = -1.5e-3, 1.5e-3
    elif spin_bh==0.: 
        ax1 = ax[1]
        tmax = 700
        #ymin,ymax = -1.e-3, 1.e-3
    elif spin_bh==0.25: 
        ax1 = ax[2]
        tmax = 900
    elif spin_bh==0.5: 
        ax1 = ax[3]
        tmax = 1000
    elif spin_bh==0.75: 
        ax1 = ax[4]
        tmax = 1200
    
    
    if (compare_nsbh[key]['type']=='NSBH'):
        datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
        tmerge_nsbh = compare_nsbh[key]['t_merger_src']
        ihspin = glob.glob(os.path.join(datadir, 'qlm_mp_j1?0?.maximum.asc'))[0]
       

    elif(compare_nsbh[key]['type']=='BBH'):
        datadir = compare_nsbh[key]['dir']
        tmerge_bbh = compare_nsbh[key]['t_merger_src'] 
        ihspin = glob.glob(os.path.join(datadir, 'qlm_mp_j1?2?.maximum.asc'))[0]
        color='b'
        
    filepath = os.path.join(datadir, 'Ylm_WEYLSCAL4::Psi4r_l2_m2_r70.00.asc')
    if not os.path.exists(filepath):
        filepath = os.path.join(datadir, 'mp_Psi4r_l2_m2_r70.00.asc')
    t22, psi4_re22, psi4_im22 = np.loadtxt(filepath, unpack=True, usecols=(0,1,2))
    amp = np.sqrt(psi4_re22**2 + psi4_im22**2)
    r=70
    
    #Read Radiated quantities data
    psi4rad_filepath = os.path.join(datadir, 'psi4radAway_l6_r70.00.asc')
        
    t_rad, Erad_rad, Kx_rad, Ky_rad, Kz_rad, Edot_rad, Pxdot_rad, Pydot_rad, Pzdot_rad, Jxdot_rad, \
    Jydot_rad, Jzdot_rad, Jx_rad, Jy_rad, Jz_rad\
             = np.loadtxt(psi4rad_filepath, unpack=True, usecols=(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14))
    
    
    Jz_adm = compare_nsbh[key]['J_adm'][-1]
    #M_adm = compare_nsbh[key]['M_adm']
    
    #Read BH Spin data
    #t_bh, sx_bh, sy_bh, sz_bh = np.loadtxt(ihspin_filepath, unpack=True, usecols=(0,1,2,3))
    R_ex_gwrad=float((os.path.basename(psi4rad_filepath).split('.asc')[0]).split('_r')[1])  
    
    t_ih, sz = np.loadtxt(ihspin, unpack=True, usecols=(1,2))
    #print("{}: ADM Angular Momentum - {:>2.3g}, Radiated Angular Momentum = {:>2.5g}, BH Spin = {:>2.5g}, \
    #BH spin + GW ang mom = {:>2.3g}".format(labelname, Jz_adm, Jz_rad[-1], sz[-1], np.absolute(Jz_rad[-1]) + np.absolute(sz[-1])))
   
    i=0

    #print(labelname,  np.absolute(J_adm)-np.absolute(Jz_rad[-1]) -np.absolute(sz[-1]))
    
    ax1.plot(t_rad, -np.absolute(Jz_rad )+ np.absolute(Jz_adm),c=color,lw=2,ls='-', label=labelname)
    ax1.plot(t_ih, sz, c='k',ls='--', label='BH Spin')
    ax1.set_xlabel(r'$T/M$')
    ax1.set_ylabel(r'$|J^{ADM} - J^{GW}|\,/M^2$', fontsize=18)
    ax1.set_xlim(0,tmax)
    ax1.legend()
    ax1.grid(True)
    loc = plticker.MultipleLocator(base=100) # this locator puts ticks at regular intervals
    
    

plt.tight_layout()
#fig.suptitle('Figure 1: BBH vs NSBH', fontsize=20, y=1.022)
#plt.savefig(os.path.join(figdir, 'Radiated_AngMom.png'), dpi=200)
plt.show()
plt.close()

i=0

ADM Energy and Angular Momentum Conservation

We next investigate the conservation of ADM energy and angular momentum to determine the energy and angular momentum trapped in the disk. Few interesting observations -

  1. While the fraction of ADM energy emitted in gravitational waves increase with BH spin, reverse trend is observed in Angular momentum.
  2. Due to stellar disruption, the mass and angular momentum of $a_BH>0.5$ continues to grow. In contrast, the remaining systems show a more BBH like behavior.
  3. The disk reaches upto 2% of ADM energy and 7% of angular momentum for highest spinning case.
In [24]:
#Energy And Angular MomentumPlot:

i=0

plt.figure(figsize=(15,10))

color = [ 'k','b','orangered','m','forestgreen']
ls = ['-', '--', ':','-.','--',':']

print('-'*100)
print('System \t\t M_bh/M_adm \t E_gw/M_adm) \t E_Disk/M_adm \t S_bh/J_adm \t J_gw/J_adm \t J_disk/J_adm')
print('-'*100)
for i, key in enumerate(sorted(compare_nsbh.keys())): #[0:3:2])):
    
    labelname = compare_nsbh[key]['label']
    col = compare_nsbh[key]['color']
    if (compare_nsbh[key]['type']=='NSBH'):
        datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
        ls='-'
        ihspin = glob.glob(os.path.join(datadir, 'qlm_mp_j1?0?.maximum.asc'))[0]
        bhdiag = os.path.join(datadir, 'BH_diagnostics.ah1.gp')
       
    elif(compare_nsbh[key]['type']=='BBH'):
        datadir = compare_nsbh[key]['dir']
        ihspin = glob.glob(os.path.join(datadir, 'qlm_mp_j1?2?.maximum.asc'))[0]
        bhdiag = os.path.join(datadir, 'BH_diagnostics.ah3.gp')
        ls='--'
        
    filepath = os.path.join(datadir, 'psi4radAway_l6_r70.00.asc')
    
    if not os.path.exists(filepath):
        print('File note found for {}'.format(labelname))
        
    t, Erad, Kx, Ky, Kz, Edot, Pxdot, Pydot, Pzdot, Jxdot, Jydot, Jzdot, Jx_rad, Jy_rad, Jz_rad\
             = np.loadtxt(filepath, unpack=True, usecols=(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14))

    
    psi4_filepath = os.path.join(datadir, 'Ylm_WEYLSCAL4::Psi4r_l2_m2_r70.00.asc')
    if not os.path.exists(psi4_filepath):
        psi4_filepath = os.path.join(datadir, 'mp_Psi4r_l2_m2_r70.00.asc')
        
    t22, psi4_re22, psi4_im22 = np.loadtxt(psi4_filepath, unpack=True, usecols=(0,1,2))
    amp = np.sqrt(psi4_re22**2 + psi4_im22**2)
    max_amp_idx = np.argmax(amp==np.amax(amp))
    
    J_adm = compare_nsbh[key]['J_adm'][-1]
    M_adm = compare_nsbh[key]['M_adm']
    
    bhdiag_filepath = os.path.join(datadir, 'BH_diagnostics.ah1.gp')
    t_ah, irr_mass = np.loadtxt(bhdiag, unpack=True, usecols=(1,26))
    t_ih, sz = np.loadtxt(ihspin, unpack=True, usecols=(1,2))
    
    spin_splinefit = splrep(t_ih[::8], sz[::8],s=0)
    sz_fit = splev(t_ah, spin_splinefit, der=0)
    
    hzn_mass = np.sqrt(irr_mass**2 + 0.25*sz_fit**2/irr_mass**2)
    
    
    print("{:15s}  {:>2.4g} \t {:>2.4g} \t {:>2.4g} \t {:>2.4g} \t {:>2.4g} \t {:>2.4g}".format(labelname,\
                                    hzn_mass[-1]/M_adm, Erad[-1]/M_adm, 1.0-(hzn_mass[-1]+ Erad[-1])/M_adm, \
                                    sz[-1]/J_adm, Jz_rad[-1]/J_adm, 1.0 - \
                                        (np.absolute(sz[-1]/J_adm) + np.absolute(Jz_rad[-1]/J_adm)) ))
    
    plt.subplot(221)
    plt.plot(t-t22[max_amp_idx], Erad/M_adm,ls=ls, c=col, label=labelname)
    plt.xlabel(r'$T/M$')
    plt.ylabel(r'$E_{rad}/M_{ADM}$')
    plt.xlim(-400,200)
    plt.legend()
    
    plt.subplot(222)
    plt.plot(t_ah-t22[max_amp_idx]+70, hzn_mass/M_adm, ls=ls, c=col, label=labelname )
    plt.xlabel(r'$T/M$')
    plt.ylabel(r'$M_{BH}/M_{ADM}$')
    plt.xlim(-400,200)
    plt.legend()
    
    plt.subplot(223)
    plt.plot(t-t22[max_amp_idx], Jz_rad/J_adm,ls=ls, c=col, label=labelname)
    plt.xlabel(r'$T/M$')
    plt.ylabel(r'${J}_{rad}/J_{ADM}$')
    plt.xlim(-150,100)
    plt.legend()
    
    plt.subplot(224)
    plt.plot(t_ih-t22[max_amp_idx], sz/J_adm,ls=ls, c=color[i], label=labelname)
    plt.xlabel(r'$T/M$')
    plt.ylabel(r'${S}_{z}/J_{ADM}$')
    plt.xlim(-150,100)
    plt.ylim(-0.5, 1)
    plt.legend()

 
print('\n')
plt.tight_layout()

#plt.savefig(os.path.join(figdir,'C3_Radiated_Energy_AngMom.png'), dpi=100)
plt.show()
plt.close()
----------------------------------------------------------------------------------------------------
System 		 M_bh/M_adm 	 E_gw/M_adm) 	 E_Disk/M_adm 	 S_bh/J_adm 	 J_gw/J_adm 	 J_disk/J_adm
----------------------------------------------------------------------------------------------------
$a_{BH} = -0.5$  0.9932 	 0.006803 	 1.997e-06 	 0.6771 	 -0.3194 	 0.00358
$a_{BH} = 0$     0.9902 	 0.009529 	 0.0002978 	 0.8234 	 -0.1706 	 0.006079
$a_{BH} = 0.25$  0.9869 	 0.01115 	 0.001902 	 0.8301 	 -0.1561 	 0.01381
$a_{BH} = 0.5$   0.9794 	 0.0121 	 0.008481 	 0.8183 	 -0.1446 	 0.03704
$a_{BH} = 0.75$  0.9664 	 0.01261 	 0.02102 	 0.7949 	 -0.1338 	 0.07131


Black Hole Kicks

Kicks vary as a function of spin but trends are not very clear. For $a_{BH}=-0.5$, black hole is strongly accelerated and is able to retain its velocity instead of dissipating in gravitational waves. This is also consistent with least energy emitted in gravitational radiation.

As the spin increases, the acceleration is quickly followed by

In [25]:
# BH Kicks
import matplotlib.gridspec as gridspec
import matplotlib as mpl
mpl.rc('lines', linewidth=2, color='r')
mpl.rc('font', size=18)
mpl.rc('axes', labelsize=18, grid=True)
mpl.rc('xtick', labelsize=14)
mpl.rc('ytick', labelsize=14)
mpl.rc('legend', fontsize=16)


color = ['k','b', 'orangered','m','coral','forestgreen']
ls = ['-', '--', ':','-.','--',':']

fig = plt.figure(constrained_layout=True, figsize=(14,15))
gs = gridspec.GridSpec(ncols=2, nrows=3, figure=fig)

ax1 = fig.add_subplot(gs[0,0])
ax2 = fig.add_subplot(gs[0,1])
ax3 = fig.add_subplot(gs[1,0])
ax4 = fig.add_subplot(gs[1,1])
ax5 = fig.add_subplot(gs[2,0])
ax6 = fig.add_subplot(gs[2,1])

for i, key in enumerate(sorted(compare_nsbh.keys())):
    
    
    labelname = compare_nsbh[key]['label']
    massratio = '5' 
    spin = float(labelname.split('=')[1][:-1])
    color = compare_nsbh[key]['color']
    if (compare_nsbh[key]['type']=='BBH'):
    
        datadir = compare_nsbh[key]['dir']
        t_merge = compare_nsbh[key]['t_merger_src']
        t_comp = t_merge + 300.
        
        bhdiag3 = os.path.join(datadir,'BH_diagnostics.ah3.gp')
        t_ah3, area_ah3,irrmass_ah3, arealrad_ah3 = np.loadtxt(bhdiag3, unpack=True, usecols=(1,25,26,27))
        
        ihspin3 = os.path.join(datadir,'ihspin_hn_2.asc')
        t_ih3, sx3, sy3, sz3, px3, py3, pz3 = np.loadtxt(ihspin3, unpack=True, usecols=(0,1,2,3,4,5,6))
        
        if spin<0.25:
            irrmass_splinefit = splrep(t_ah3[::8],irrmass_ah[::8],s=0)
            spinz_splinefit = splrep(t_ih3[::2], sz[::2], s=0)
        else :
            irrmass_splinefit = splrep(t_ah3[::128],irrmass_ah[::128],s=0)
            spinz_splinefit = splrep(t_ih3[::16], sz[::16], s=0)
    
        
        sz_fit = splev(t_ah3, spinz_splinefit, der=0)
        hzn_mass = np.sqrt(irrmass_ah3**2 + 0.25*sz_fit**2/irrmass_ah3**2)
        hzn_mass_fit = splrep(t_ah3, hzn_mass, s=0)
        
        
        psi4rad_filepath = os.path.join(datadir, 'psi4radAway_l6_r70.00.asc')
   
        trad, Erad, Kxrad, Kyrad, Kzrad, Edotrad, Pxdotrad, Pydotrad, Pzdotrad\
             = np.loadtxt(psi4rad_filepath, unpack=True, usecols=(0,1,2,3,4,5,6,7,8))
    
        hzn_mass_evalfit = splev(trad, hzn_mass_fit, der=0)
        
        ax1.plot(trad - (t_merge+70), Kxrad, c=color[i], ls='--')
        ax1.plot(t_ih3, px, c=color[i], ls=':')
        ax2.semilogy(trad - t_merge-70, Pxdotrad, c=color[i], ls='--')

        ax3.plot(trad - (t_merge+70), Kyrad, c=color[i], ls='--')
        ax3.plot(t_ih3, py, c=color[i], ls=':')
        ax4.semilogy(trad - t_merge-70, Pydotrad, c=color[i], ls='--')

        ax5.plot(trad - (t_merge+70), np.sqrt(Kxrad**2 + Kyrad**2 + Kzrad**2), c=color[i], ls='--')
        ax5.plot(t_ih3, np.sqrt(px**2 + py**2 + pz**2), c=color[i], ls=':')
       
        ax6.semilogy(trad - t_merge-70, np.sqrt(Pxdotrad**2 + Pydotrad**2 + Pzdotrad**2), c=color[i], ls='--')

        print(labelname,bbh_irrmass_bh3, bbh_spinmag_bh3[-1])
        
    elif(compare_nsbh[key]['type']=='NSBH'):
        
        datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
        t_merge = compare_nsbh[key]['t_merger_src']
        t_comp = t_merge + 300.
        
        
        bhdiag = os.path.join(datadir,'BH_diagnostics.ah1.gp')
        t_ah, area_ah,irrmass_ah, arealrad_ah = np.loadtxt(bhdiag, unpack=True, usecols=(1,25,26,27))
        
        ihspin = os.path.join(datadir,'ihspin_hn_0.asc')
        t_ih, sx, sy, sz, px, py, pz = np.loadtxt(ihspin, unpack=True, usecols=(0,1,2,3,4,5,6))
     
        
        if spin<0.25:
            irrmass_splinefit = splrep(t_ah[::8],irrmass_ah[::8],s=0)
            spinz_splinefit = splrep(t_ih[::2], sz[::2], s=0)
        else :
            irrmass_splinefit = splrep(t_ah[::128],irrmass_ah[::128],s=0)
            spinz_splinefit = splrep(t_ih[::16], sz[::16], s=0)
    
        
        sz_fit = splev(t_ah, spinz_splinefit, der=0)
        hzn_mass = np.sqrt(irrmass_ah**2 + 0.25*sz_fit**2/irrmass_ah**2)
        hzn_mass_fit = splrep(t_ah, hzn_mass, s=0)
        
        
        psi4rad_filepath = os.path.join(datadir, 'psi4radAway_l6_r70.00.asc')
    
        trad, Erad, Kxrad, Kyrad, Kzrad, Edotrad, Pxdotrad, Pydotrad, Pzdotrad\
             = np.loadtxt(psi4rad_filepath, unpack=True, usecols=(0,1,2,3,4,5,6,7,8))
        hzn_mass_evalfit = splev(trad-70, hzn_mass_fit, der=0)
        
        ax1.plot(trad - (t_merge+70), Kxrad/hzn_mass_evalfit, c=color, ls='-',label=labelname)
        #ax1.plot(t_ih - t_merge, px, c=color, ls=':')
        ax2.plot(trad - t_merge-70, Pxdotrad, c=color, ls='-',label=labelname)

        ax3.plot(trad - (t_merge+70), Kyrad/hzn_mass_evalfit, c=color, ls='-',label=labelname)
        #ax3.plot(t_ih - t_merge, py, c=color, ls=':')
        ax4.plot(trad - t_merge-70, Pydotrad, c=color, ls='-',label=labelname)

        ax5.plot(trad - (t_merge+70), np.sqrt(Kxrad**2 + Kyrad**2 + Kzrad**2)/hzn_mass_evalfit, c=color, ls='-',label=labelname)
        #ax5.plot(t_ih - t_merge, np.sqrt(px**2 + py**2 + pz**2), c=color, ls=':')
       
        ax6.semilogy(trad - t_merge-70, np.sqrt(Pxdotrad**2 + Pydotrad**2 + Pzdotrad**2), c=color, ls='-',label=labelname)

        
        
        
ax1.set_xlabel(r'$T/M$', fontsize=20)
ax1.set_ylabel(r'$ K_x$', fontsize=20)
ax2.set_xlabel(r'$T/M$', fontsize=20)
ax2.set_ylabel(r'$\dot{K}_x$', fontsize=20)
ax3.set_xlabel(r'$T/M$', fontsize=20)
ax3.set_ylabel(r'$K_y$', fontsize=20)
ax4.set_xlabel(r'$T/M$', fontsize=20)
ax4.set_ylabel(r'$\dot{K}_y$', fontsize=20)
ax5.set_xlabel(r'$T/M$', fontsize=20)
ax5.set_ylabel(r'$K/M$', fontsize=20)
ax6.set_xlabel(r'$T/M$', fontsize=20)
ax6.set_ylabel(r'$\dot{K}$', fontsize=20)

ax1.set_xlim(-50,80)
ax2.set_xlim(-50,80)
ax3.set_xlim(-50,80)
ax4.set_xlim(-50,80)
ax5.set_xlim(-50,80)
ax6.set_xlim(-50,80)

#ax2.set_ylim(2e-8,2e-4)
#ax4.set_ylim(2e-8,2e-4)
ax6.set_ylim(2e-8,2e-4)

ax1.legend(fontsize=12, loc=4)
ax2.legend(fontsize=12, loc=1)
ax3.legend(fontsize=12, loc=4)
ax4.legend(fontsize=12)
ax5.legend(fontsize=12)
ax6.legend(fontsize=12)
    
ax1.grid(True)
ax2.grid(True)
ax3.grid(True)
ax4.grid(True)
ax5.grid(True)
ax6.grid(True)

plt.tight_layout()
#plt.savefig(os.path.join(figdir, 'BH_Kick.png'), dpi=100)
plt.show()
plt.close()
This figure was using constrained_layout==True, but that is incompatible with subplots_adjust and or tight_layout: setting constrained_layout==False. 

QuasiLocal Quantities

Black Hole and Neutron Star Mass

In [26]:
#BH Horizon Mass Plots - Comparison of mass  and spin of final black hole at 300M after merger
color = ['r', 'b', 'k','m','coral','forestgreen']
ls = ['-', '--', ':','-.','--',':']

bbh_irrmass_dict = {}
bbh_hznmass_dict = {}

print("System \t \t Time \t Init. Irr Mass    Final Irr Mass \t Final Spin(a) \tFinal Hzn Mass   \n ")
print('-'*115)


for i, key in enumerate(sorted(compare_nsbh.keys())):
    
    
    labelname = compare_nsbh[key]['label']
    massratio = '5' 
    spin = float(labelname.split('=')[1][:-1])
    
    if (compare_nsbh[key]['type']=='BBH'):
        datadir = compare_nsbh[key]['dir']
        t_merge = compare_nsbh[key]['t_merger_src']
        t_comp = t_merge + 300.
        
        bhdiag1 = os.path.join(datadir,'BH_diagnostics.ah1.gp')
        t_ah1, area_ah1,irrmass_ah1, arealrad_ah1 = np.loadtxt(bhdiag1, unpack=True, usecols=(1,25,26,27))
        
        ihspin1 = os.path.join(datadir,'ihspin_hn_0.asc')
        t_ih1, sx1, sy1, sz1 = np.loadtxt(ihspin1, unpack=True, usecols=(0,1,2,3))
        
        bbh_irrmass_bh1 = irrmass_ah1[0]
        bbh_spinmag_bh1 = (sx1**2 + sy1**2  + sz1**2)**0.5
        bbh_hznmass_bh1 = (bbh_irrmass_bh1**2 + 0.25*bbh_spinmag_bh1[0]**2/bbh_irrmass_bh1**2)**0.5
        
        
        bhdiag3 = os.path.join(datadir,'BH_diagnostics.ah3.gp')
        t_ah3, area_ah3,irrmass_ah3, arealrad_ah3 = np.loadtxt(bhdiag3, unpack=True, usecols=(1,25,26,27))
        
        ihspin3 = os.path.join(datadir,'ihspin_hn_2.asc')
        t_ih3, sx3, sy3, sz3 = np.loadtxt(ihspin3, unpack=True, usecols=(0,1,2,3))
        
        #Set all the simulations to same final time before comparing mass and spins
        area_ah3       = area_ah3    [t_ah3 < t_comp]
        irrmass_ah3    = irrmass_ah3 [t_ah3 < t_comp]
        arealrad_ah3   = arealrad_ah3[t_ah3 < t_comp]
        t_ah3          = t_ah3       [t_ah3 < t_comp]
        
        sx3   = sx3  [t_ih3 < t_comp]
        sy3   = sy3  [t_ih3 < t_comp]
        sz3   = sz3  [t_ih3 < t_comp]
        t_ih3 = t_ih3[t_ih3 < t_comp]
        
        
        bbh_irrmass_bh3 = irrmass_ah3[-1]
        bbh_spinmag_bh3 = (sx3**2 + sy3**2  + sz3**2)**0.5
        bbh_hznmass_bh3 = (bbh_irrmass_bh3**2 + 0.25*bbh_spinmag_bh3[-1]**2/bbh_irrmass_bh3**2)**0.5
        bbh_dimspin_bh3 = bbh_spinmag_bh3[-1]/bbh_hznmass_bh3**2
        
        bbh_irrmass_dict[massratio]  = bbh_irrmass_bh3
        bbh_hznmass_dict[massratio]  = bbh_hznmass_bh3
        print("%-10s \t %g \t  %g \t\t %g \t\t %g \t %g \t %g "%(labelname,bbh_irrmass_bh1, bbh_irrmass_bh3, \
                                            bbh_irrmass_bh3 - bbh_irrmass_bh3, bbh_dimspin_bh3,\
                                            bbh_hznmass_bh3,  bbh_hznmass_bh3 - bbh_hznmass_bh3))
        
        
    elif(compare_nsbh[key]['type']=='NSBH'):
       
        datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
        t_merge = compare_nsbh[key]['t_merger_src']
        t_comp = t_merge + 400.
        
        
        bhdiag = os.path.join(datadir,'BH_diagnostics.ah1.gp')
        t_ah, area_ah,irrmass_ah, arealrad_ah = np.loadtxt(bhdiag, unpack=True, usecols=(1,25,26,27))
        
        ihspin = os.path.join(datadir,'ihspin_hn_0.asc')
        t_ih, sx, sy, sz = np.loadtxt(ihspin, unpack=True, usecols=(0,1,2,3))
     
        mass_growth_rate = np.gradient(irrmass_ah[::8])/np.gradient(t_ah[::8])
        
        
        #Set all the simulations to same final time before comparing mass and spins
        area_ah       = area_ah    [t_ah<t_comp]
        irrmass_ah    = irrmass_ah [t_ah<t_comp]
        arealrad_ah   = arealrad_ah[t_ah<t_comp]
        t_ah          = t_ah       [t_ah<t_comp]
       
        sx   = sx  [t_ih<t_comp]
        sy   = sy  [t_ih<t_comp]
        sz   = sz  [t_ih<t_comp]
        t_ih = t_ih[t_ih<t_comp]
        
        
        
        nsbh_irrmass_bh1 = irrmass_ah[0]
        nsbh_spinmag_bh1 = (sx**2 + sy**2  + sz**2)**0.5
        nsbh_hznmass_bh1 = (nsbh_irrmass_bh1**2 + 0.25*nsbh_spinmag_bh1[0]**2/nsbh_irrmass_bh1**2)**0.5
        
        nsbh_irrmass_bh3 = irrmass_ah[-1]
        nsbh_spinmag_bh3 = (sx**2 + sy**2  + sz**2)**0.5
        nsbh_hznmass_bh3 = (nsbh_irrmass_bh3**2 + 0.25*nsbh_spinmag_bh3[-1]**2/nsbh_irrmass_bh3**2)**0.5
        nsbh_dimspin_bh3 = nsbh_spinmag_bh3[-1]/nsbh_hznmass_bh3**2
            
        print("%-1s \t %.1f \t %g \t  %g \t \t%g \t %g "%(labelname[1:-1], t_ah[-1]-t_merge, nsbh_irrmass_bh1,nsbh_irrmass_bh3, \
                                                  nsbh_dimspin_bh3, nsbh_hznmass_bh3))
       
        
System 	 	 Time 	 Init. Irr Mass    Final Irr Mass 	 Final Spin(a) 	Final Hzn Mass   
 
-------------------------------------------------------------------------------------------------------------------
a_{BH} = -0.5 	 344.4 	 0.804443 	  0.98568 	 	0.119759 	 0.987458 
a_{BH} = 0 	 400.0 	 0.833313 	  0.960979 	 	0.421161 	 0.984137 
a_{BH} = 0.25 	 400.0 	 0.826648 	  0.936397 	 	0.567263 	 0.980658 
a_{BH} = 0.5 	 400.0 	 0.804443 	  0.898644 	 	0.708301 	 0.973027 
a_{BH} = 0.75 	 400.0 	 0.755167 	  0.837836 	 	0.845233 	 0.956544 
In [42]:
#Masss andd Spin Plot (Paper)
import matplotlib.gridspec as gridspec
import matplotlib as mpl
mpl.rc('lines', linewidth=2, color='r')
mpl.rc('font', size=18)
mpl.rc('axes', labelsize=18, grid=True)
mpl.rc('xtick', labelsize=14)
mpl.rc('ytick', labelsize=14)
mpl.rc('legend', fontsize=16)


color = ['k','b', 'orangered','m','coral','forestgreen']
ls = ['-', '--', ':','-.','--',':']

fig = plt.figure(constrained_layout=True, figsize=(14,10))
gs = gridspec.GridSpec(ncols=2, nrows=2, figure=fig)

ax1 = fig.add_subplot(gs[0,0])
ax2 = fig.add_subplot(gs[0,1])
ax3 = fig.add_subplot(gs[1,0])
ax4 = fig.add_subplot(gs[1,1])
#ax5 = fig.add_subplot(gs[2,0])
#ax6 = fig.add_subplot(gs[2,1])

for i, key in enumerate(sorted(compare_nsbh.keys())):
    
    
    labelname = compare_nsbh[key]['label']
    massratio = '5' 
    spin = float(labelname.split('=')[1][:-1])
    color = compare_nsbh[key]['color']
    if (compare_nsbh[key]['type']=='BBH'):
    
        datadir = compare_nsbh[key]['dir']
        t_merge = compare_nsbh[key]['t_merger_src']
        t_comp = t_merge + 300.
        
        bhdiag3 = os.path.join(datadir,'BH_diagnostics.ah3.gp')
        t_ah3, area_ah3,irrmass_ah3, arealrad_ah3 = np.loadtxt(bhdiag3, unpack=True, usecols=(1,25,26,27))
        
        ihspin3 = os.path.join(datadir,'ihspin_hn_2.asc')
        t_ih3, sx3, sy3, sz3,px3,py3,pz3 = np.loadtxt(ihspin3, unpack=True, usecols=(0,1,2,3,4,5,6))
        
        #Set all the simulations to same final time before comparing mass and spins
        irrmass_ah3    = irrmass_ah3 [t_ah3 < t_comp]
        t_ah3          = t_ah3       [t_ah3 < t_comp]
        
        sx3   = sx3  [t_ih3 < t_comp]
        sy3   = sy3  [t_ih3 < t_comp]
        sz3   = sz3  [t_ih3 < t_comp]
        t_ih3 = t_ih3[t_ih3 < t_comp]
        
        
        bbh_irrmass_bh3 = irrmass_ah3[-1]
        bbh_spinmag_bh3 = (sx3**2 + sy3**2  + sz3**2)**0.5
        bbh_hznmass_bh3 = (bbh_irrmass_bh3**2 + 0.25*bbh_spinmag_bh3[-1]**2/bbh_irrmass_bh3**2)**0.5
        bbh_dimless_spin_bh3 = bbh_spinmag_bh3[-1]/bbh_hznmass_bh3**2
        
        ax1.axhline(y=bbh_irrmass_bh3,c=color[i], ls='--')
        ax3.axhline(y=bbh_spinmag_bh3[-1],c=color[i], ls='--')
       
        psi4rad_filepath = os.path.join(datadir, 'psi4radAway_l6_r70.00.asc')
   
        trad, Erad, Kxrad, Kyrad, Kzrad, Edotrad, Pxdotrad, Pydotrad, Pzdotrad\
             = np.loadtxt(psi4rad_filepath, unpack=True, usecols=(0,1,2,3,4,5,6,7,8))
        #ax5.plot(trad - (t_merge+70), np.sqrt(Kxrad**2 + Kyrad**2 + Kzrad**2), c=color[i], ls='--')
        #ax5.plot(t_ih3, np.sqrt(px**2 + py**2 + pz**2), c=color[i], ls=':')
       
        #ax6.semilogy(trad - t_merge-70, np.sqrt(Pxdotrad**2 + Pydotrad**2 + Pzdotrad**2), c=color[i], ls='--')

        print(labelname,bbh_irrmass_bh3, bbh_spinmag_bh3[-1])
    elif(compare_nsbh[key]['type']=='NSBH'):
        
        datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
        t_merge = compare_nsbh[key]['t_merger_src']
        t_comp = t_merge + 300.
        
        
        bhdiag = os.path.join(datadir,'BH_diagnostics.ah1.gp')
        t_ah, area_ah,irrmass_ah, arealrad_ah = np.loadtxt(bhdiag, unpack=True, usecols=(1,25,26,27))
        
        ihspin = os.path.join(datadir,'ihspin_hn_0.asc')
        t_ih, sx, sy, sz = np.loadtxt(ihspin, unpack=True, usecols=(0,1,2,3))
     
        
        if spin<0.25:
            irrmass_splinefit = splrep(t_ah[::8],irrmass_ah[::8],s=0)
            spinz_splinefit = splrep(t_ih[::2], sz[::2], s=0)
        else :
            irrmass_splinefit = splrep(t_ah[::128],irrmass_ah[::128],s=0)
            spinz_splinefit = splrep(t_ih[::16], sz[::16], s=0)
    
        #Compute mass growth rate on every 8th iteration to reduce the high frequency noise in the derivatives 
        #Truncation and round off errors leads to high frequency noise which gets amplified on computing the derivatives
        #Computing the derivative on a larger time step averages out these effect to zero)
       
        mass_growth_rate_fit =  splev(t_ah, irrmass_splinefit, der=1) 
        mass_growth_rate = np.gradient(irrmass_ah[::32])/np.gradient(t_ah[::32])
        
        spin_growth_rate_fit =  splev(t_ih, spinz_splinefit, der=1) 
        spin_growth_rate = np.gradient(sz)/np.gradient(t_ih)
        
        sz_fit = splev(t_ah, spinz_splinefit, der=0)
        hzn_mass = np.sqrt(irrmass_ah**2 + 0.25*sz_fit**2/irrmass_ah**2)
        hzn_mass_fit = splrep(t_ah, hzn_mass, s=0)
        
        
        ax1.plot(t_ah-t_merge, irrmass_ah-irrmass_ah[0], c=color, ls='-',label=labelname)
        #ax1.plot(t_ah-t_merge, splev(t_ah, irrmass_splinefit), color[i], ls='--')
        
        ax2.semilogy((t_ah-t_merge), mass_growth_rate_fit, c=color, ls='-',label=labelname)
        #ax2.semilogy((t_ah-t_merge)[::32], mass_growth_rate, c=color[i], ls='--',label=labelname[:2])
        
        ax3.plot(t_ih-t_merge, sz-sz[0], c=color, ls='-',label=labelname)
        ax4.semilogy(t_ih-t_merge, spin_growth_rate, c=color, ls='-',lw=2,label=labelname)
        #ax4.semilogy(t_ih-t_merge, spin_growth_rate_fit, c=color[i-3], ls='--', lw=2,label=labelname[:2])
        
        psi4rad_filepath = os.path.join(datadir, 'psi4radAway_l6_r70.00.asc')
    
        trad, Erad, Kxrad, Kyrad, Kzrad, Edotrad, Pxdotrad, Pydotrad, Pzdotrad\
             = np.loadtxt(psi4rad_filepath, unpack=True, usecols=(0,1,2,3,4,5,6,7,8))
        hzn_mass_evalfit = splev(trad, hzn_mass_fit, der=0)
        
        #ax5.plot(trad-(t_merge+70), np.sqrt(Kxrad**2 + Kyrad**2 + Kzrad**2)/hzn_mass_evalfit, c=color, ls='-',label=labelname)
        #ax5.semilogy(trad, np.sqrt(Px**2 + Py**2 + Pz**2), c=color[i], ls='--',label=labelname[:2])

        #ax6.semilogy(trad-t_merge-70, np.sqrt(Pxdotrad**2 + Pydotrad**2 + Pzdotrad**2), c=color, ls='-',label=labelname)
        
        
        
ax1.set_xlabel(r'$T/M$', fontsize=20)
ax1.set_ylabel(r'$ \Delta M_i/M$', fontsize=20)
ax2.set_xlabel(r'$T/M$', fontsize=20)
ax2.set_ylabel(r'$\dot{M}_i$', fontsize=20)
ax3.set_xlabel(r'$T/M$', fontsize=20)
ax3.set_ylabel(r'$\Delta S_z/M^2$', fontsize=20)
ax4.set_xlabel(r'$T/M$', fontsize=20)
ax4.set_ylabel(r'$\dot{S}_z/M$', fontsize=20)
#ax5.set_xlabel(r'$T/M$', fontsize=20)
#ax5.set_ylabel(r'$|P/M|$', fontsize=20)
#ax6.set_xlabel(r'$T/M$', fontsize=20)
#ax6.set_ylabel(r'$|\dot{P}|$', fontsize=20)

ax1.set_xlim(-200,300)
ax2.set_xlim(-50,80)
ax3.set_xlim(-200,300)
ax4.set_xlim(-50,80)
#ax5.set_xlim(-200,300)
#ax6.set_xlim(-50,80)

#ax1.set_ylim(-0.01,0.005)
ax2.set_ylim(2e-5,0.1)
ax4.set_ylim(2e-5,0.2)
#ax6.set_ylim(2e-8,2e-4)

ax1.legend(loc=4)
ax2.legend()
ax3.legend()
ax4.legend()
#ax5.legend()
#ax6.legend()

#plt.axhline(y=irrmass_ah[0], ls='--', c='k')
#plt.ylim(0.65, 0.92)
#ax5.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
    
ax1.grid(True)
ax2.grid(True)
ax3.grid(True)
ax4.grid(True)
#ax5.grid(True)
#ax6.grid(True)
plt.tight_layout()
#plt.savefig(os.path.join(figdir, 'BH_MassPlots.png'), dpi=100)
plt.show()
plt.close()

Rest Mass outside BH

In [31]:
#Rest Mass outside the horizon
#subplt_num=1
print("System \t\t Time - T_merger\t   Mass (R>1.5)/M_* ")
print('-'*100)
        
plt.figure(figsize=(8,8))
for i, key in enumerate(sorted(compare_nsbh.keys())):
    
    labelname = compare_nsbh[key]['label']
    massratio = '5' 
    spin = float(labelname.split('=')[1][:-1])
    t_merge = compare_nsbh[key]['t_merger_src']
    color = compare_nsbh[key]['color']
    if (compare_nsbh[key]['type']=='BBH'):
        continue
    
    elif(compare_nsbh[key]['type']=='NSBH'):
       
        datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
        
        #Volume Integrals Data
        vi_grmhd = os.path.join(datadir,'volume_integrals-GRMHD.asc')
        time, m_ns, m_bh3, m_0d7, m_1d4, m_20, m_35  = np.loadtxt(vi_grmhd, unpack=True, usecols=(0,5,6,7,8,9,10))
        
        
        #Choose merger point as reference for calculation
        time_centered = time-t_merge
        idx = np.where(time_centered<400)  
        
        time_centered = time_centered[idx]
        
        if compare_nsbh[key]['symmetry']:
            m_ns  = 2.*m_ns[idx] #accounting for symmetries
            m_bh3 = 2.*m_bh3[idx] #accounting for symmetries
            m_0d7 = 2.*m_0d7[idx] #accounting for symmetries
            m_1d4 = 2.*m_1d4[idx] #accounting for symmetries
            m_20  = 2.*m_20[idx] #accounting for symmetries
            m_35  = 2.*m_35[idx] #accounting for symmetries
            
            
        #This is approximate mass and we are ignoring the mass in sphere of 0.7M due to lack of data. This mass
        #will be of order of 1e-15 and hence can be neglected
        
        restmass_atm = m_0d7[0] - m_ns[0]
        #print('Atmosphere rest mass = %g (M) = %g (M_*), m_ns=%g'%(restmass_atm, restmass_atm/m_ns[0], m_ns[0]))
        
           
        mass_outside_hzn = m_1d4[idx] - restmass_atm
        mass_outside_R35 = m_35[idx]  - m_35[0]
        mass_outside_R20 = m_20[idx]  - m_20[0]
        m_disk_R35_M =  mass_outside_hzn - mass_outside_R35
        m_disk_R20_M =  mass_outside_hzn - mass_outside_R20
        m_disk_R35_Mstar = m_disk_R35_M/m_ns[0]
        m_disk_R20_Mstar = m_disk_R20_M/m_ns[0]
        m_ejecta_R35_M = mass_outside_R35
        mass_outhzn_Mstar =  mass_outside_hzn/m_ns[0]
        
        print("%-10s \t %.2f \t \t %.2g "%(labelname[1:-1], time_centered[-1], mass_outhzn_Mstar[-1]))
        
          

        plt.semilogy(time_centered, mass_outhzn_Mstar, c=color, ls='-', label=labelname)  
        #plt.semilogy(time_centered, m_disk_R35_Mstar, c=color[i-3], ls='--', label=labelname[:2]+' (R = 35M)')
        
        plt.xlabel(r'$T/M$')
        plt.ylabel(r'$M_{td}/M_o$')
        plt.ylim(1e-5,2)
        plt.xlim(-100,+400)
        #plt.tick_params(left=True, right=True, which='both')
        #plt.axhline(y=0.06)
        plt.legend()
        
        
plt.grid(True)
#plt.savefig(os.path.join(figdir, 'ReviewFigures/RevFig_RestMassOutsideBH.png'), dpi=100)
plt.show()
plt.close()
System 		 Time - T_merger	   Mass (R>1.5)/M_* 
----------------------------------------------------------------------------------------------------
a_{BH} = -0.5 	 344.39 	 	 7.8e-05 
a_{BH} = 0 	 400.00 	 	 0.00064 
a_{BH} = 0.25 	 400.00 	 	 0.0052 
a_{BH} = 0.5 	 400.00 	 	 0.035 
a_{BH} = 0.75 	 400.00 	 	 0.13 
In [ ]:
#Foucart Fits. -

# arXiv:1807.00011 - Equation 4

alpha_fit = 0.406
beta_fit = 0.139
gamma_fit = 0.255
delta_fit = 1.761


def remnant_mass_fit(nsbh_dict):
    
    
    comp = nsbh_dict['comp']
    q = nsbh_dict['q']
    a_bh = nsbh_dict['a_bh']
    eta = q/(1. + q)**2.
    
    z1 = 1. + ((1. - a_bh**2)**(1./3))*((1. + a_bh)**(1./3) + (1. - a_bh)**(1./3))
    z2 = np.sqrt(3.*a_bh**2. + z1**2.)
    
    isco_radius_hat = 3. + z2 - np.sign(a_bh)*np.sqrt((3. - z1)*(3. + z1 + 2*z2))
    
    remnant_mass = (max(alpha_fit*(1-2.0*comp)/(eta**(1./3)) - beta_fit*isco_radius_hat*comp/eta + gamma_fit,\
                        0))**delta_fit
    return remnant_mass


q2_dict = {'comp':0.1529, 'q':2, 'a_bh':0.} #0.683}
q3_dict = {'comp':0.1536, 'q':3, 'a_bh':0.} #0.563}
q5_dict = {'comp':0.1543, 'q':5, 'a_bh':0.} #0.4203}
q7_dict = {'comp':0.145 , 'q':5, 'a_bh':0.0}

print('Q2: Remnant mass = %g'%remnant_mass_fit(q2_dict))
print('Q3: Remnant mass = %g'%remnant_mass_fit(q3_dict))
print('Q5: Remnant mass = %g'%remnant_mass_fit(q5_dict))
print('Q7 (Test): Remnant mass = %g'%remnant_mass_fit(q7_dict))

Multipole Moments

Here, we look at the multipole moments. The $M^i_{mp}$ refers to mass multipoles while $J^i_{mp}$ refers to spin multipoles. For mass multipoles, the odd moments are 0 while for spin multipoles, the even moments are zero. The zeroeth mass multipole is basically the Christodoulou mass, as can be seen in the first plot and zeroeth spin multipole in the spin angular momentum of black hole.

Things to look into

  • Comparison between BBH and NSBH mass multipoles.
  • Does accreting matter affect the multipole moments of black hole?
  • Why do higher moments continue to exist (i.e. nonzero) for final black hole?
  • Read - Isolated and dynamical horizons and their applications, Section 5.3.
In [46]:
#Quasi Local Moments
color = ['r', 'b', 'k','m','coral','forestgreen']
ls = ['-', '--', ':','-.','--',':']

bbh_ah3_mass = 0

#print("System \t\t\t Initial BH Mass \t Final BH Mass \t NSBH - BBH \n ")

fig1, ax1 = plt.subplots(2,2,figsize=(15,8))
ax1 = np.concatenate(ax1)

fig2, ax2 = plt.subplots(2,2,figsize=(15,8))
ax2 = np.concatenate(ax2)

for i, key in enumerate(sorted(compare_nsbh.keys())):
   
    labelname = compare_nsbh[key]['label']
    spin = float(labelname.split('=')[1][:-1])
    t_merge = compare_nsbh[key]['t_merger_src']
    color = compare_nsbh[key]['color']
    if (compare_nsbh[key]['type']=='BBH'):
        datadir = compare_nsbh[key]['dir']
        continue
        
    elif(compare_nsbh[key]['type']=='NSBH'):
        datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
        

    qlm_mass_files = sorted(glob.glob(os.path.join(datadir, 'qlm_mp_m*.norm2.asc')))
    if len(qlm_mass_files)>0:
        bhdiag = os.path.join(datadir,'BH_diagnostics.ah1.gp')
        t_ah, area_ah,irrmass_ah, arealrad_ah = np.loadtxt(bhdiag, unpack=True, usecols=(1,25,26,27))
        t_ah = np.around(t_ah, decimals=4)
        
        ihspin = os.path.join(datadir,'ihspin_hn_0.asc')
        t_ih, sx,sy,sz = np.loadtxt(ihspin, unpack=True, usecols=(0,1,2,3))
        t_ih = np.around(t_ih, decimals=4)
        
        t_common = np.intersect1d(t_ah, t_ih)
        ahf_idx = np.isin(t_ah,t_common)
        ih_idx = np.isin(t_ih,t_common)
        
        spinmag = np.sqrt(sx**2 + sy**2  + sz**2)
        hznmass_ah = (irrmass_ah[ahf_idx]**2 + 0.25*spinmag[ih_idx]**2/irrmass_ah[ahf_idx]**2)**0.5
        
        
    for qlm_file in qlm_mass_files:
        
        mass_mp_label = int((os.path.basename(qlm_file).split('[')[0])[-1])
        t, mass_moment = np.loadtxt(qlm_file, unpack=True, usecols=(1,2))
        
        if mass_mp_label%2==1 or mass_mp_label>7: continue
            
        ax = ax1[mass_mp_label//2]
        
        ax.plot(t-t_merge,mass_moment, c=color, label='%s'%(labelname))
        ax.set_xlabel('Time')
        ax.set_ylabel(r'$M_{}$'.format(mass_mp_label))
        ax.legend()
        ax.grid(True)
        #ax.set_ylim(-0.4,0.4)
        ax.set_xlim(-50,60)
        
    
    #if len(qlm_mass_files)>0: ax11.plot(t_ah[ahf_idx], hznmass_ah, ls='--', c='k',label='%s:BH Mass'%(labelname))
    
    
    qlm_am_files = sorted(glob.glob(os.path.join(datadir, 'qlm_mp_j*.norm2.asc')))
    
    for qlm_file in qlm_am_files:
        
        am_mp_label = int((os.path.basename(qlm_file).split('[')[0])[-1])
        if am_mp_label%2==0: continue
        t, j_moment = np.loadtxt(qlm_file, unpack=True, usecols=(1,2))
        ax =  ax2[(am_mp_label-1)//2]
        ax.plot(t-t_merge, j_moment, c=color,label='%s'%(labelname))
        ax.set_xlabel('Time')
        ax.set_ylabel(r'$J_{}$'.format(am_mp_label))
        ax.legend()
        ax.grid(True)
        ax.set_ylim(-0.3,1)
        ax.set_xlim(-50,70)
    

fig1.tight_layout()
fig1.show()
#plt.close()

fig2.tight_layout()
fig2.show()
#plt.close()
                        

Other resources -

Additional Plots

In [ ]:
simdir1 = '/nethome/numrel/datafiles/Waveforms/HR-series/q5.00/D10_q5.00_a-0.71_0.00_m240/'


psi4_filepath = os.path.join(simdir1, 'mp_WeylScal4::Psi4i_l2_m2_r70.00.asc')
t_psi4, re_psi4, im_psi4 = np.loadtxt(psi4_filepath, unpack=True, usecols=(0,1,2)) 
amp = np.sqrt(re_psi4**2 + im_psi4**2)

t_max_amp = t_psi4[amp == np.amax(amp)]
psi4rad_filepath = os.path.join(simdir1, 'psi4analysis_r70.00.asc') #'psi4radAway_l6_r70.00.asc')
trad, Erad, Kxrad, Kyrad, Kzrad, Edotrad, Pxdotrad, Pydotrad, Pzdotrad\
             = np.loadtxt(psi4rad_filepath, unpack=True, usecols=(0,1,2,3,4,5,6,7,8))
             
Kmag = np.sqrt(Kxrad**2 + Kyrad**2 + Kzrad**2)
plt.plot(trad - t_max_amp, Kmag)


simdir2 = '/nethome/numrel/datafiles/Waveforms/HR-series/q5.00/D10_q5.00_a0.6_0.0_m240/'
psi4_filepath = os.path.join(simdir2, 'mp_WeylScal4::Psi4i_l2_m2_r70.00.asc')
t_psi4, re_psi4, im_psi4 = np.loadtxt(psi4_filepath, unpack=True, usecols=(0,1,2)) 
amp = np.sqrt(re_psi4**2 + im_psi4**2)

t_max_amp = t_psi4[amp == np.amax(amp)]

psi4rad_filepath = os.path.join(simdir2, 'psi4analysis_r70.00.asc') #'psi4radAway_l6_r70.00.asc')
trad, Erad, Kxrad, Kyrad, Kzrad, Edotrad, Pxdotrad, Pydotrad, Pzdotrad\
             = np.loadtxt(psi4rad_filepath, unpack=True, usecols=(0,1,2,3,4,5,6,7,8))
             
Kmag = np.sqrt(Kxrad**2 + Kyrad**2 + Kzrad**2)
plt.plot(trad - t_max_amp, Kmag)
plt.xlim(-50,100)
plt.show()

Orbital Frequency

In [ ]:
# Orbital Frequency and rate of change of separation - 
plt.figure(figsize=(14,8))
print("\nNote - Mixed binary derivatives are extremely noise. Hence, to tackle the noise, we use spline interpolation \
to fit the original data and compute the derivatives using these fitting functions. While they do provide a \
good approximation to original fits (Tested with BBH), they fail to capture the high frequency oscillations in rdot and omega \
for mixed binaries")
print('-------------'*9)
print('\n\n')
for i, key in enumerate(sorted(compare_nsbh.keys())):
    
    t = orbdyn_dict[key]['time']
   
    label = compare_nsbh[key]['label']
            
    #plt.plot(orbdyn_dict[key]['time'], orbdyn_dict[key]['rdot'], ls=compare_nsbh[key]['ls'],c=compare_nsbh[key]['color'],label=label)
    plt.plot(orbdyn_dict[key]['time'], splev(orbdyn_dict[key]['time'], orbdyn_dict[key]['sepfit'], der=1), ls=compare_nsbh[key]['ls'], c=compare_nsbh[key]['color'],label=label+'Spline Fit')
    
    plt.xlabel('Time', fontsize=16)
    plt.ylabel('dr/dt', fontsize=16)
    plt.legend()
plt.title(r'Figure 3 - $dR_{sep}/dt$', fontsize=20)
plt.xlim(0,900)
plt.ylim(-0.025,0.015)
plt.show()
plt.close()

print('\n\n')
    
fig = plt.figure(figsize=(15,7))

for i,key in enumerate(sorted(compare_nsbh.keys())):
    
    label = compare_nsbh[key]['label']

    omega_fit = splev(orbdyn_dict[key]['time'], orbdyn_dict[key]['phasefit'], der=1)
    
    if 'Q2-NSBH' in label:
        time = orbdyn_dict[key]['time'][:-140]
        omega_fit = omega_fit[:-140]
        
    else:
        time = orbdyn_dict[key]['time']

    #print("%s: Initial Orbital Frequency = %g"%(label, orbdyn_dict[key]['omega'][0]))
    plt.subplot(121)

    #plt.plot(orbdyn_dict[key]['time'],orbdyn_dict[key]['omega'],ls=compare_nsbh[key]['ls'],c=compare_nsbh[key]['color'],label=label)
    plt.plot(time, omega_fit, ls=compare_nsbh[key]['ls'],c=compare_nsbh[key]['color'],label=label+' Spline Derivative')

    plt.xlabel('Time/M', fontsize=16)
    plt.ylabel(r'$M\Omega_{orb}$', fontsize=16)
    plt.ylim(0.022,0.4)
    plt.legend()
    
    
    plt.subplot(122)
        
    #plt.plot(orbdyn_dict[key]['time'], orbdyn_dict[key]['omega']/2./np.pi, ls=compare_nsbh[key]['ls'],c=compare_nsbh[key]['color'],label=label)
    plt.plot(time, omega_fit, ls=compare_nsbh[key]['ls'],c=compare_nsbh[key]['color'],label=label+' Spline Derivative')
    
    plt.xlabel('Time/M', fontsize=16)
    plt.ylabel(r'$M\Omega_{orb}$ ', fontsize=16)
    plt.legend()
    plt.xlim(0,600)
    #plt.ylim(0.022,0.1)
    
plt.tight_layout()
fig.suptitle('Figure 4 - Orbital Frequency', y=1.08, fontsize=20)
#plt.show()
plt.close()

Gravitational Waves

In [ ]:
#Re(Psi4) and Amplitude plots: Thesis Figures
import matplotlib.ticker as plticker
from matplotlib.legend_handler import HandlerLine2D, HandlerTuple



i=0

nsbh_color = ['k','b','r','darkviolet','cadetblue', 'r','k','m','darkolivegreen','coral','forestgreen', 'blueviolet']
bbh_color = ['dimgray','royalblue','orangered','darkorchid', 'darkcyan', 'teal','orangered', 'grey',]
ls = ['--', '-', ':','-.','--',':']

fig, (ax1,ax2) = plt.subplots(2, 1, figsize=(10,8))
 

for i, key in enumerate(sorted(compare_nsbh.keys())):
    labelname = compare_nsbh[key]['label']
    
    if 'Q2' not in labelname: continue
    if 'Q2' in labelname: 
        q=2
        tmax = 900
        l = 0.242
        b = 0.892
        w = 0.26
        h = 0.065
    if 'Q3' in labelname: 
        q=3
        tmax = 1000
        l = 0.228
        b = 0.89
        w = 0.235
        h = 0.067
        
    if 'Q5' in labelname: 
        q=5
        tmax = 1200
        l = 0.203
        w = 0.213
        h = 0.0706
        b = 0.224
    
    
    ls = compare_nsbh[key]['ls']
    if (compare_nsbh[key]['type']=='NSBH'):
        datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
        tmerge_nsbh = compare_nsbh[key]['t_merger_src']
        color=nsbh_color
        ls='-'
        lw=1.5
        ihspin_filepath = os.path.join(datadir, 'ihspin_hn_0.asc')
        
    elif(compare_nsbh[key]['type']=='BBH'):
        datadir = compare_nsbh[key]['dir']
        tmerge_bbh = compare_nsbh[key]['t_merger_src']
        color=bbh_color 
        ls='--'
        lw=2
        ihspin_filepath = os.path.join(datadir, 'ihspin_hn_2.asc')
        
        
    #Read Psi4 data
    psi4_filepath = os.path.join(datadir, 'Ylm_WEYLSCAL4::Psi4r_l2_m2_r70.00.asc')
    if not os.path.exists(filepath):
        psi4_filepath = os.path.join(datadir, 'mp_Psi4r_l2_m2_r70.00.asc')
        
    t22, psi4_re22, psi4_im22 = np.loadtxt(psi4_filepath, unpack=True, usecols=(0,1,2))
    amp = np.sqrt(psi4_re22**2 + psi4_im22**2)
    R_ex_gw=float((os.path.basename(psi4_filepath).split('.asc')[0]).split('_r')[1])  
    
    
    #Read Radiated quantities data
    psi4rad_filepath = os.path.join(datadir, 'psi4radAway_l6_r70.00.asc')
        
    t_rad, Erad_rad, Kx_rad, Ky_rad, Kz_rad, Edot_rad, Pxdot_rad, Pydot_rad, Pzdot_rad, Jxdot_rad, Jydot_rad, Jzdot_rad, Jx_rad, Jy_rad, Jz_rad\
             = np.loadtxt(psi4rad_filepath, unpack=True, usecols=(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14))
    
    #Read BH Spin data
    t_bh, sx_bh, sy_bh, sz_bh = np.loadtxt(ihspin_filepath, unpack=True, usecols=(0,1,2,3))
    R_ex_gwrad=float((os.path.basename(psi4_filepath).split('.asc')[0]).split('_r')[1])  
    
    i=0

    ax1.plot(t22, r*psi4_re22,c=color[i],ls=ls, lw=lw)#, label=labelname)
    labelplot = ax1.plot([],[],c='k', ls=ls, label=labelname[3:])
    #ax1.plot(t22,amp,ls=ls,c=color)
    #ax1.set_xlabel(r'$T/M$')
    ax1.set_ylabel(r'$Re(\psi_4)\,r\,M$')
    ax1.set_xlim(0,tmax)
    ax1.set_ylim(-6e-2,6e-2)
    ax1.legend(loc=3)
    ax1.grid(True)
    ax1.axvline(150, ls='--', c='grey')
    ax1.axvline(451, ls='--', c='grey')

    ax1.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
    loc = plticker.MultipleLocator(base=100) # this locator puts ticks at regular intervals
    ax1.xaxis.set_major_locator(loc)
    
    i+=1
    ax2.plot(t_rad, Jz_rad,c=color[i],ls=ls, lw=lw)
    if 'NSBH' in labelname:
        labelplot1 = ax2.plot([],[],c=color[i], ls=ls,  label=r'$J_z$')
        labelplot2 = ax2.plot([],[],c=color[i+1], ls=ls, label=r'$S_z$')
    
    ax2.set_xlabel(r'$T/M$')
    if 'NSBH' in labelname: ax2.set_ylabel(r'$J_z/M^2$')#,color=color[i])
    ax2.set_xlim(0,tmax)
    ax2.grid(True)
    ax2.legend(loc=3)
    
    #ax2.axvline(150, ls='--', c='grey')
    #ax2.axvline(450, ls='--', c='grey')
    ax2.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
    loc = plticker.MultipleLocator(base=100) # this locator puts ticks at regular intervals
    ax2.xaxis.set_major_locator(loc)
    
    #Twin plot in second plot
    i+=1
    ax21 = ax2.twinx()  # instantiate a second axes that shares the same x-axis
    
    if 'NSBH' in labelname:
        Sz_plot = ax21.plot(t_bh, sz_bh, lw=lw, color=color[i], label=r'$S_z$')
        #ax21.yaxis.set_ticklabels([])
        ax21.set_ylabel(r'$S_z/M^2$')#, color=color[i])  # we already handled the x-label with ax1
    
        #ax21.legend()#, numpoints=1,handler_map={tuple: HandlerTuple(ndivide=None)})
    
    else:
        print(sz_bh[-1])
        Sz_plot = ax21.axhline(y = sz_bh[-1], c=color[i], lw=lw, ls=ls)
        #ax21.yaxis.set_ticklabels([])
        ax21.tick_params(axis='y', labelright='off',right=False)
        ax21.minorticks_off()
       
    ax21.set_ylim(-0.05,0.7)
    #ax21.yaxis.set_ticklabels(np.around(np.arange(0, 0.62,0.15), decimals=2))
    #    
    loc = plticker.MultipleLocator(base=0.15) # this locator puts ticks at regular intervals
    ax21.yaxis.set_major_locator(loc)
       
    plt.grid(False)
    
    
    
    #Inset in Plot 1
    i=0
    left, bottom, width, height = [l, b, w, h]   
    ax12 = fig.add_axes([left, bottom, width, height])

    min_idx = np.amin(np.argwhere(t22>=150))
    max_idx = np.amin(np.argwhere(t22>=451))
    
    ax12.plot(t22, psi4_re22*r,c=color[i],ls=ls, lw=lw)
    ax12.axhline(y=0,c='grey', alpha=0.5, lw=1)
    ax12.axvline(x=200,c='grey', alpha=0.5, lw=1)
    ax12.axvline(x=300,c='grey', alpha=0.5, lw=1)
    ax12.axvline(x=400,c='grey', alpha=0.5, lw=1)
    ax12.set_xlim(150, 450)
    ax12.set_ylim(-1.3e-3, 1.3e-3)
    ax12.yaxis.tick_right()
    #ax12.tick_params(axis='both',labelbottom='off', labelright='off',bottom=False,right=False)
    ax12.grid(b=True, which='major', axis='x')
    ax12.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
    ax12.tick_params(axis='y', labelsize=11)
    ax12.yaxis.offsetText.set_fontsize(11)
    
    ax12.minorticks_off()
    ax2.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
    plt.xticks([])
    #plt.yticks([])
    

plt.tight_layout()
#fig.suptitle('Figure 1: BBH vs NSBH', fontsize=20, y=1.022)
print(os.path.join(figdir, 'Psi4Plots_q%d.png'%q))
plt.savefig(os.path.join(figdir, 'Psi4Plots_q%d.png'%q), dpi=200)
plt.show()
plt.close()
    
In [ ]:
#Re(Psi4) and Amplitude plots: Thesis Figures
import matplotlib.ticker as plticker
from matplotlib.legend_handler import HandlerLine2D, HandlerTuple



i=0

nsbh_color = ['k','b','r','darkviolet','cadetblue', 'r','k','m','darkolivegreen','coral','forestgreen', 'blueviolet']
bbh_color = ['dimgray','royalblue','orangered','darkorchid', 'darkcyan', 'teal','orangered', 'grey',]
ls = ['--', '-', ':','-.','--',':']

fig, (ax1,ax2) = plt.subplots(2, 1, figsize=(10,8))
 

for i, key in enumerate(sorted(compare_nsbh.keys())):
    labelname = compare_nsbh[key]['label']
    
    if 'Q3' not in labelname: continue
    if 'Q2' in labelname: 
        q=2
        tmax = 900
        l = 0.242
        b = 0.883
        w = 0.26
        h = 0.065
    if 'Q3' in labelname: 
        q=3
        tmax = 1000
        l = 0.228
        b = 0.89
        w = 0.235
        h = 0.067
        
    if 'Q5' in labelname: 
        q=5
        tmax = 1200
        l = 0.203
        w = 0.213
        h = 0.0706
        b = 0.224
    
    
    ls = compare_nsbh[key]['ls']
    if (compare_nsbh[key]['type']=='NSBH'):
        datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
        tmerge_nsbh = compare_nsbh[key]['t_merger_src']
        color=nsbh_color
        ls='-'
        lw=1.5
        ihspin_filepath = os.path.join(datadir, 'ihspin_hn_0.asc')
        
    elif(compare_nsbh[key]['type']=='BBH'):
        datadir = compare_nsbh[key]['dir']
        tmerge_bbh = compare_nsbh[key]['t_merger_src']
        color=bbh_color 
        ls='--'
        lw=2
        ihspin_filepath = os.path.join(datadir, 'ihspin_hn_2.asc')
        
        
    #Read Psi4 data
    psi4_filepath = os.path.join(datadir, 'Ylm_WEYLSCAL4::Psi4r_l2_m2_r70.00.asc')
    if not os.path.exists(filepath):
        psi4_filepath = os.path.join(datadir, 'mp_Psi4r_l2_m2_r70.00.asc')
        
    t22, psi4_re22, psi4_im22 = np.loadtxt(psi4_filepath, unpack=True, usecols=(0,1,2))
    amp = np.sqrt(psi4_re22**2 + psi4_im22**2)
    R_ex_gw=float((os.path.basename(psi4_filepath).split('.asc')[0]).split('_r')[1])  
    
    
    #Read Radiated quantities data
    psi4rad_filepath = os.path.join(datadir, 'psi4radAway_l6_r70.00.asc')
        
    t_rad, Erad_rad, Kx_rad, Ky_rad, Kz_rad, Edot_rad, Pxdot_rad, Pydot_rad, Pzdot_rad, Jxdot_rad, Jydot_rad, Jzdot_rad, Jx_rad, Jy_rad, Jz_rad\
             = np.loadtxt(psi4rad_filepath, unpack=True, usecols=(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14))
    
    #Read BH Spin data
    t_bh, sx_bh, sy_bh, sz_bh = np.loadtxt(ihspin_filepath, unpack=True, usecols=(0,1,2,3))
    R_ex_gwrad=float((os.path.basename(psi4_filepath).split('.asc')[0]).split('_r')[1])  
    
    i=0

    ax1.plot(t22, r*psi4_re22,c=color[i],ls=ls, lw=lw)#, label=labelname)
    labelplot = ax1.plot([],[],c='k', ls=ls, label=labelname[3:])
    #ax1.plot(t22,amp,ls=ls,c=color)
    #ax1.set_xlabel(r'$T/M$')
    ax1.set_ylabel(r'$Re(\psi_4)\,r\,M$')
    ax1.set_xlim(0,tmax)
    ax1.legend(loc=3)
    ax1.grid(True)
    ax1.axvline(150, ls='--', c='grey')
    ax1.axvline(451, ls='--', c='grey')
    ax1.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
    loc = plticker.MultipleLocator(base=100) # this locator puts ticks at regular intervals
    ax1.xaxis.set_major_locator(loc)
    
    i+=1
    ax2.plot(t_rad, Jz_rad,c=color[i],ls=ls, lw=lw)
    if 'NSBH' in labelname:
        labelplot1 = ax2.plot([],[],c=color[i], ls=ls,  label=r'$J_z$')
        labelplot2 = ax2.plot([],[],c=color[i+1], ls=ls, label=r'$S_z$')
    
    ax2.set_xlabel(r'$T/M$')
    if 'NSBH' in labelname: ax2.set_ylabel(r'$J_z/M^2$')#,color=color[i])
    ax2.set_xlim(0,tmax)
    ax2.grid(True)
    ax2.legend(loc=3)
    
    #ax2.axvline(150, ls='--', c='grey')
    #ax2.axvline(450, ls='--', c='grey')
    ax2.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
    loc = plticker.MultipleLocator(base=100) # this locator puts ticks at regular intervals
    ax2.xaxis.set_major_locator(loc)
    
    #Twin plot in second plot
    i+=1
    ax21 = ax2.twinx()  # instantiate a second axes that shares the same x-axis
    
    if 'NSBH' in labelname:
        Sz_plot = ax21.plot(t_bh, sz_bh, lw=lw, color=color[i], label=r'$S_z$')
        #ax21.yaxis.set_ticklabels([])
        ax21.set_ylabel(r'$S_z/M^2$')#, color=color[i])  # we already handled the x-label with ax1
    
        #ax21.legend()#, numpoints=1,handler_map={tuple: HandlerTuple(ndivide=None)})
    
    else:
        print(sz_bh[-1])
        Sz_plot = ax21.axhline(y = sz_bh[-1], c=color[i], lw=lw, ls=ls)
        #ax21.yaxis.set_ticklabels([])
        ax21.tick_params(axis='y', labelright='off',right=False)
        ax21.minorticks_off()
       
    ax21.set_ylim(-0.05,0.6)
    #ax21.yaxis.set_ticklabels(np.around(np.arange(0, 0.62,0.15), decimals=2))
    #    
    loc = plticker.MultipleLocator(base=0.15) # this locator puts ticks at regular intervals
    ax21.yaxis.set_major_locator(loc)
       
    plt.grid(False)
    
    
    
    #Inset in Plot 1
    i=0
    left, bottom, width, height = [l, b, w, h]   
    ax12 = fig.add_axes([left, bottom, width, height])

    min_idx = np.amin(np.argwhere(t22>=150))
    max_idx = np.amin(np.argwhere(t22>=451))
    
    ax12.plot(t22, psi4_re22*r,c=color[i],ls=ls, lw=lw)
    ax12.axhline(y=0,c='grey', alpha=0.5, lw=1)
    ax12.axvline(x=200,c='grey', alpha=0.5, lw=1)
    ax12.axvline(x=300,c='grey', alpha=0.5, lw=1)
    ax12.axvline(x=400,c='grey', alpha=0.5, lw=1)
    ax12.set_xlim(150, 450)
    ax12.set_ylim(-1e-3, 1e-3)
    ax12.yaxis.tick_right()
    #ax12.tick_params(axis='both',labelbottom='off', labelright='off',bottom=False,right=False)
    ax12.grid(b=True, which='major', axis='x')
    ax12.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
    ax12.tick_params(axis='y', labelsize=11)
    ax12.yaxis.offsetText.set_fontsize(11)
    
    ax12.minorticks_off()
    
    ax2.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
    plt.xticks([])
    #plt.yticks([])
    

plt.tight_layout()
#fig.suptitle('Figure 1: BBH vs NSBH', fontsize=20, y=1.022)
print(os.path.join(figdir, 'Psi4Plots_q%d.png'%q))
plt.savefig(os.path.join(figdir, 'Psi4Plots_q%d.png'%q), dpi=200)
plt.show()
plt.close()
    
In [ ]:
#Re(Psi4) and Amplitude plots: Thesis Figures
import matplotlib.ticker as plticker
from matplotlib.legend_handler import HandlerLine2D, HandlerTuple



i=0

nsbh_color = ['k','mediumblue','r','darkviolet','cadetblue', 'r','k','m','darkolivegreen','coral','forestgreen', 'blueviolet']
bbh_color = ['dimgray','royalblue','tomato','darkorchid', 'darkcyan', 'teal','orangered', 'grey',]
ls = ['--', '-', ':','-.','--',':']

fig, (ax1,ax2) = plt.subplots(2, 1, figsize=(10,8))
 

for i, key in enumerate(sorted(compare_nsbh.keys())):
    labelname = compare_nsbh[key]['label']
    
    if 'Q5' not in labelname: continue
    if 'Q2' in labelname: 
        q=2
        tmax = 900
        l = 0.242
        b = 0.883
        w = 0.26
        h = 0.065
    if 'Q3' in labelname: 
        q=3
        tmax = 1000
        l = 0.228
        b = 0.89
        w = 0.235
        h = 0.067
        
    if 'Q5' in labelname: 
        q=5
        tmax = 1200
        l = 0.2105
        b = 0.8916
        w = 0.1983
        h = 0.065
    
    
    ls = compare_nsbh[key]['ls']
    if (compare_nsbh[key]['type']=='NSBH'):
        datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
        tmerge_nsbh = compare_nsbh[key]['t_merger_src']
        color=nsbh_color
        ls='-'
        lw=1.5
        ihspin_filepath = os.path.join(datadir, 'ihspin_hn_0.asc')
        
    elif(compare_nsbh[key]['type']=='BBH'):
        datadir = compare_nsbh[key]['dir']
        tmerge_bbh = compare_nsbh[key]['t_merger_src']
        color=bbh_color 
        ls='--'
        lw=2
        ihspin_filepath = os.path.join(datadir, 'ihspin_hn_2.asc')
        
        
    #Read Psi4 data
    psi4_filepath = os.path.join(datadir, 'Ylm_WEYLSCAL4::Psi4r_l2_m2_r70.00.asc')
    if not os.path.exists(filepath):
        psi4_filepath = os.path.join(datadir, 'mp_Psi4r_l2_m2_r70.00.asc')
        
    t22, psi4_re22, psi4_im22 = np.loadtxt(psi4_filepath, unpack=True, usecols=(0,1,2))
    amp = np.sqrt(psi4_re22**2 + psi4_im22**2)
    R_ex_gw=float((os.path.basename(psi4_filepath).split('.asc')[0]).split('_r')[1])  
    
    
    #Read Radiated quantities data
    psi4rad_filepath = os.path.join(datadir, 'psi4radAway_l6_r70.00.asc')
        
    t_rad, Erad_rad, Kx_rad, Ky_rad, Kz_rad, Edot_rad, Pxdot_rad, Pydot_rad, Pzdot_rad, Jxdot_rad, Jydot_rad, Jzdot_rad, Jx_rad, Jy_rad, Jz_rad\
             = np.loadtxt(psi4rad_filepath, unpack=True, usecols=(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14))
    
    #Read BH Spin data
    t_bh, sx_bh, sy_bh, sz_bh = np.loadtxt(ihspin_filepath, unpack=True, usecols=(0,1,2,3))
    R_ex_gwrad=float((os.path.basename(psi4_filepath).split('.asc')[0]).split('_r')[1])  
    
    i=0

    ax1.plot(t22, r*psi4_re22,c=color[i],ls=ls, lw=lw)#, label=labelname)
    labelplot = ax1.plot([],[],c='k', ls=ls, label=labelname[3:])
    #ax1.plot(t22,amp,ls=ls,c=color)
    #ax1.set_xlabel(r'$T/M$')
    ax1.set_ylabel(r'$Re(\psi_4)\,r\,M$')
    ax1.set_xlim(0,tmax)
    ax1.set_ylim(-3e-2, 3e-2)
    ax1.legend(loc=3)
    ax1.grid(True)
    ax1.axvline(150, ls='--', c='grey')
    ax1.axvline(450, ls='--', c='grey')
    ax1.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
    loc = plticker.MultipleLocator(base=100) # this locator puts ticks at regular intervals
    ax1.xaxis.set_major_locator(loc)
    
    i+=1
    ax2.plot(t_rad, Jz_rad,c=color[i],ls=ls, lw=lw)
    if 'NSBH' in labelname:
        labelplot1 = ax2.plot([],[],c=color[i], ls=ls,  label=r'$J_z$')
        labelplot2 = ax2.plot([],[],c=color[i+1], ls=ls, label=r'$S_z$')
    
    ax2.set_xlabel(r'$T/M$')
    if 'NSBH' in labelname: ax2.set_ylabel(r'$J_z/M^2$')#,color=color[i])
    ax2.set_xlim(0,tmax)
    ax2.grid(True)
    ax2.legend(loc=3)
    
    #ax2.axvline(150, ls='--', c='grey')
    #ax2.axvline(450, ls='--', c='grey')
    ax2.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
    loc = plticker.MultipleLocator(base=100) # this locator puts ticks at regular intervals
    ax2.xaxis.set_major_locator(loc)
    
    #Twin plot in second plot
    i+=1
    ax21 = ax2.twinx()  # instantiate a second axes that shares the same x-axis
    
    if 'NSBH' in labelname:
        Sz_plot = ax21.plot(t_bh, sz_bh, lw=lw, color=color[i], label=r'$S_z$')
        #ax21.yaxis.set_ticklabels([])
        ax21.set_ylabel(r'$S_z/M^2$')#, color=color[i])  # we already handled the x-label with ax1
    
        #ax21.legend()#, numpoints=1,handler_map={tuple: HandlerTuple(ndivide=None)})
    
    else:
        print(sz_bh[-1])
        Sz_plot = ax21.axhline(y = sz_bh[-1], c=color[i], lw=lw, ls=ls)
        #ax21.yaxis.set_ticklabels([])
        ax21.tick_params(axis='y', labelright='off',right=False)
        ax21.minorticks_off()
       
    ax21.set_ylim(-0.05,0.5)
    #ax21.yaxis.set_ticklabels(np.around(np.arange(0, 0.62,0.15), decimals=2))
    #    
    loc = plticker.MultipleLocator(base=0.1) # this locator puts ticks at regular intervals
    ax21.yaxis.set_major_locator(loc)
       
    plt.grid(False)
    
    
    
    #Inset in Plot 1
    i=0
    left, bottom, width, height = [l, b, w, h]   
    ax12 = fig.add_axes([left, bottom, width, height])

    min_idx = np.amin(np.argwhere(t22>=150))
    max_idx = np.amin(np.argwhere(t22>=451))
    
    ax12.plot(t22, psi4_re22*r,c=color[i],ls=ls, lw=lw)
    ax12.axhline(y=0,c='grey', alpha=0.5, lw=1)
    ax12.axvline(x=200,c='grey', alpha=0.5, lw=1)
    ax12.axvline(x=300,c='grey', alpha=0.5, lw=1)
    ax12.axvline(x=400,c='grey', alpha=0.5, lw=1)
    ax12.set_xlim(150, 450)
    ax12.set_ylim(-1e-3, 1e-3)
    ax12.yaxis.tick_right()
    #ax12.tick_params(axis='both',labelbottom='off', labelright='off',bottom=False,right=False)
    ax12.grid(b=True, which='major', axis='x')
    ax12.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
    ax12.tick_params(axis='y', labelsize=11)
    ax12.yaxis.offsetText.set_fontsize(11)
    
    ax12.minorticks_off()
    
    ax2.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
    plt.xticks([])
    #plt.yticks([])
    

plt.tight_layout()
#fig.suptitle('Figure 1: BBH vs NSBH', fontsize=20, y=1.022)
print(os.path.join(figdir, 'Psi4Plots_q%d.png'%q))
plt.savefig(os.path.join(figdir, 'Psi4Plots_q%d.png'%q), dpi=200)
plt.show()
plt.close()
    
In [ ]:
# Q2 psi4 frequency oscillations
from scipy import signal

plt.figure(figsize=(15,6))
for i, key in enumerate(sorted(compare_nsbh.keys())):
    
    labelname = compare_nsbh[key]['label']
    if (compare_nsbh[key]['type']=='NSBH'):
        datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
        tmerge_nsbh = compare_nsbh[key]['t_merger_src']
    elif(compare_nsbh[key]['type']=='BBH'):
        datadir = compare_nsbh[key]['dir']
        continue
    
    filepath = os.path.join(datadir, 'Ylm_WEYLSCAL4::Psi4r_l2_m2_r70.00.asc')
    if not os.path.exists(filepath):
        filepath = os.path.join(datadir, 'mp_Psi4r_l2_m2_r70.00.asc')
    
    q = float((labelname.split('Q')[1]).split('-')[0])
    mtotal = 1.35*(1+q)
   
    
    if 'Q2' not in labelname:continue
    t22, psi4_re22, psi4_im22 = np.loadtxt(filepath, unpack=True, usecols=(0,1,2))

    #Compute GW frequency
    phi = np.unwrap(np.arctan2(psi4_im22,psi4_re22))
    amp = np.sqrt(psi4_re22**2 + psi4_im22**2)
    max_amp_idx = np.argmax(amp==np.amax(amp))
    gw_freq = np.gradient(phi)/np.gradient(t22)/2./np.pi
    
    #Limit the data to a small time interval during early inspiral say 200-350M
    t_low = 200
    t_high = 400
    idx_fft = np.intersect1d(np.where(t22>t_low), np.where(t22<t_high))
    time_fft = t22[idx_fft]
    
    #Fit over this data
    gw_freq_fit = np.poly1d(np.polyfit(t22[idx_fft], gw_freq[idx_fft], deg=2))
    tfreq_fit = np.arange(t_low, t_high, t22[1]-t22[0])
    
    gw_osc = (gw_freq[idx_fft]-gw_freq_fit(time_fft))/gw_freq_fit(time_fft)
    
    
    #Apply high pass filter to remove noise 
    sos = signal.butter(10, 6000*(mtotal/s1_in_Msun), 'lp',  output='sos', )
    filtered = signal.sosfilt(sos, gw_osc)
    
    
    
    plt.subplot(121)
    plt.plot(t22, gw_freq,ls='-', c=color[i], label=labelname)
    plt.plot(tfreq_fit, gw_freq_fit(tfreq_fit),ls='--', c='k', label=labelname)
    plt.xlabel(r'$Time - T_{max}$')
    plt.ylabel('Psi4 Frequency Oscillations')
    plt.xlim(100,620)
    #plt.ylim(-0.005,0.005)
    plt.ylim(-0.015,-0.01)
    
    
    #plt.subplot(222)
    #plt.plot(time_fft, gw_osc,ls=ls[i], c=color[i], label=labelname)
    #plt.plot(time_fft, filtered,ls='--', c='k', label=labelname)
    #plt.xlabel(r'$Time - T_{max}$')
    #plt.ylabel('Psi4 Frequency Oscillations')
    #plt.xlim(250,350)
    ##plt.ylim(-0.015,-0.01)
    #plt.ylim(-0.1,0.1)
    
   
    #Remove the first and last 20M of the fit to avoid side effects of polynomial fitting
    idx_fft = np.intersect1d(np.where(t22>t_low+10), np.where(t22<t_high-30))
    time_fft = t22[idx_fft]
    gw_osc = (gw_freq[idx_fft]-gw_freq_fit(time_fft))/gw_freq_fit(time_fft)
    
    #Apply high pass filter to remove noise 
    sos = signal.butter(10, 6000*(mtotal/s1_in_Msun), 'lp',  output='sos', )
    filtered = signal.sosfilt(sos, gw_osc)
   
    osc_fft = np.fft.fft(filtered)
    fft_freq = np.fft.fftfreq(np.size(time_fft), d=time_fft[2]-time_fft[1])
    fft_freq_hz = fft_freq/(mtotal/s1_in_Msun)

   
    
    
    plt.subplot(122)
    plt.plot(fft_freq_hz, np.abs(osc_fft.real+1.j*osc_fft.imag),ls=ls[i], c=color[i], label=labelname)
    plt.xlabel(r'$FFT Freq$')
    plt.ylabel('FFT(Psi4 Frequency)')
    plt.axvline(x=1880)
    plt.xlim(0,3000)
    #plt.ylim(-3,2.5)
    
plt.tight_layout()
plt.legend()
plt.show()
plt.close()
In [ ]:
# Q3 GW frequency oscillations
from scipy import signal

plt.figure(figsize=(15,6))
for i, key in enumerate(sorted(compare_nsbh.keys())):
    
    labelname = compare_nsbh[key]['label']
    if (compare_nsbh[key]['type']=='NSBH'):
        datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
        tmerge_nsbh = compare_nsbh[key]['t_merger_src']
    elif(compare_nsbh[key]['type']=='BBH'):
        datadir = compare_nsbh[key]['dir']
        continue
    
    filepath = os.path.join(datadir, 'Ylm_WEYLSCAL4::Psi4r_l2_m2_r70.00.asc')
    if not os.path.exists(filepath):
        filepath = os.path.join(datadir, 'mp_Psi4r_l2_m2_r70.00.asc')
    
    q = float((labelname.split('Q')[1]).split('-')[0])
    mtotal = 1.35*(1+q)
   
    
    if 'Q3' not in labelname:continue
    t22, psi4_re22, psi4_im22 = np.loadtxt(filepath, unpack=True, usecols=(0,1,2))

    #Compute GW frequency
    phi = np.unwrap(np.arctan2(psi4_im22,psi4_re22))
    amp = np.sqrt(psi4_re22**2 + psi4_im22**2)
    max_amp_idx = np.argmax(amp==np.amax(amp))
    gw_freq = np.gradient(phi)/np.gradient(t22)/2./np.pi
    
    #Limit the data to a small time interval during early inspiral say 200-300M
    t_low = 200
    t_high = 400
    idx_fft = np.intersect1d(np.where(t22>t_low), np.where(t22<t_high))
    time_fft = t22[idx_fft]
    
    #Fit over this data
    gw_freq_fit = np.poly1d(np.polyfit(t22[idx_fft], gw_freq[idx_fft], deg=2))
    tfreq_fit = np.arange(t_low, t_high, t22[1]-t22[0])
    
    gw_osc = (gw_freq[idx_fft]-gw_freq_fit(time_fft))/gw_freq_fit(time_fft)
    
    
    #Apply high pass filter to remove noise 
    sos = signal.butter(10, 6000*(mtotal/s1_in_Msun), 'lp',  output='sos', )
    filtered = signal.sosfilt(sos, gw_osc)
    
    
    
    plt.subplot(121) #221)
    plt.plot(t22, gw_freq,ls='-', c=color[i], label=labelname)
    plt.plot(tfreq_fit, gw_freq_fit(tfreq_fit),ls='--', c='k', label=labelname)
    plt.xlabel(r'$Time - T_{max}$')
    plt.ylabel('Psi4 Frequency Oscillations')
    plt.xlim(200,400)
    #plt.ylim(-0.005,0.005)
    plt.ylim(-0.015,-0.01)
    
    
    #plt.subplot(222)
    #plt.plot(time_fft, gw_osc,ls=ls[i], c=color[i], label=labelname)
    #plt.plot(time_fft, filtered,ls='--', c='k', label=labelname)
    #plt.xlabel(r'$Time - T_{max}$')
    #plt.ylabel('Psi4 Frequency Oscillations')
    #plt.xlim(220,380)
    ##plt.ylim(-0.015,-0.01)
    #plt.ylim(-0.1,0.1)
    
   
    #Remove the first and last 20M of the fit to avoid side effects of polynomial fitting
    idx_fft = np.intersect1d(np.where(t22>t_low+30), np.where(t22<t_high-30))
    time_fft = t22[idx_fft]
    gw_osc = (gw_freq[idx_fft]-gw_freq_fit(time_fft))/gw_freq_fit(time_fft)
    
    #Apply high pass filter to remove noise 
    sos = signal.butter(10, 6000*(mtotal/s1_in_Msun), 'lp',  output='sos', )
    filtered = signal.sosfilt(sos, gw_osc)
   
    osc_fft = np.fft.fft(filtered)
    fft_freq = np.fft.fftfreq(np.size(time_fft), d=time_fft[2]-time_fft[1])
    fft_freq_hz = fft_freq/(mtotal/s1_in_Msun)

   
    
    
    plt.subplot(122) #223)
    plt.plot(fft_freq_hz, np.abs(osc_fft.real),ls=ls[i], c=color[i], label=labelname)
    plt.xlabel(r'$FFT Freq$')
    plt.ylabel('FFT(Psi4 Frequency)')
    #plt.axvline(x=0.028)
    plt.xlim(500,2000)
    #plt.ylim(-3,2.5)
    
plt.tight_layout()
plt.legend()
plt.show()
plt.close()
In [ ]:
# Q5 GW frequency oscillations
from scipy import signal

plt.figure(figsize=(15,6))
for i, key in enumerate(sorted(compare_nsbh.keys())):
    
    labelname = compare_nsbh[key]['label']
    if (compare_nsbh[key]['type']=='NSBH'):
        datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
        tmerge_nsbh = compare_nsbh[key]['t_merger_src']
    elif(compare_nsbh[key]['type']=='BBH'):
        datadir = compare_nsbh[key]['dir']
        continue
    
    filepath = os.path.join(datadir, 'Ylm_WEYLSCAL4::Psi4r_l2_m2_r70.00.asc')
    if not os.path.exists(filepath):
        filepath = os.path.join(datadir, 'mp_Psi4r_l2_m2_r70.00.asc')
    
    q = float((labelname.split('Q')[1]).split('-')[0])
    mtotal = 1.35*(1+q)
   
    
    if 'Q5' not in labelname:continue
    t22, psi4_re22, psi4_im22 = np.loadtxt(filepath, unpack=True, usecols=(0,1,2))

    #Compute GW frequency
    phi = np.unwrap(np.arctan2(psi4_im22,psi4_re22))
    amp = np.sqrt(psi4_re22**2 + psi4_im22**2)
    max_amp_idx = np.argmax(amp==np.amax(amp))
    gw_freq = np.gradient(phi)/np.gradient(t22)/2./np.pi
    
    #Limit the data to a small time interval during early inspiral say 200-300M
    t_low = 200
    t_high = 400
    idx_fft = np.intersect1d(np.where(t22>t_low), np.where(t22<t_high))
    time_fft = t22[idx_fft]
    
    #Fit over this data
    gw_freq_fit = np.poly1d(np.polyfit(t22[idx_fft], gw_freq[idx_fft], deg=2))
    tfreq_fit = np.arange(t_low, t_high, t22[1]-t22[0])
    
    gw_osc = (gw_freq[idx_fft]-gw_freq_fit(time_fft))/gw_freq_fit(time_fft)
    
    
    #Apply high pass filter to remove noise 
    sos = signal.butter(10, 6000*(mtotal/s1_in_Msun), 'lp',  output='sos', )
    filtered = signal.sosfilt(sos, gw_osc)
    
    
    
    plt.subplot(121)   #   221)
    plt.plot(t22, gw_freq,ls='-', c=color[i], label=labelname)
    plt.plot(tfreq_fit, gw_freq_fit(tfreq_fit),ls='--', c='k', label=labelname)
    plt.xlabel(r'$Time - T_{max}$')
    plt.ylabel('Psi4 Frequency Oscillations')
    plt.xlim(t_low, t_high)
    #plt.ylim(-0.005,0.005)
    plt.ylim(-0.015,-0.01)
    
    
    #plt.subplot(222)
    #plt.plot(time_fft, gw_osc,ls=ls[i], c=color[i], label=labelname)
    #plt.plot(time_fft, filtered,ls='--', c='k', label=labelname)
    #plt.xlabel(r'$Time - T_{max}$')
    #plt.ylabel('Psi4 Frequency Oscillations')
    #plt.axvline(x=481)
    #plt.axvline(x=511)
    #plt.xlim(t_low+20,t_high-20)
    ##plt.ylim(-0.015,-0.01)
    #plt.ylim(-0.1,0.1)
    
   
    #Remove the first and last 20M of the fit to avoid side effects of polynomial fitting
    idx_fft = np.intersect1d(np.where(t22>t_low+30), np.where(t22<t_high-30))
    time_fft = t22[idx_fft]
    gw_osc = (gw_freq[idx_fft]-gw_freq_fit(time_fft))/gw_freq_fit(time_fft)
    
    #Apply high pass filter to remove noise 
    sos = signal.butter(10, 6000*(mtotal/s1_in_Msun), 'lp',  output='sos', )
    filtered = signal.sosfilt(sos, gw_osc)
   
    osc_fft = np.fft.fft(filtered)
    fft_freq = np.fft.fftfreq(np.size(time_fft), d=time_fft[2]-time_fft[1])
    fft_freq_hz = fft_freq/(mtotal/s1_in_Msun)

   
    
    
    plt.subplot(122)        #223)
    plt.plot(fft_freq_hz, np.abs(osc_fft.real),ls=ls[i], c=color[i], label=labelname)
    plt.xlabel(r'$FFT Freq$')
    plt.ylabel('FFT(Psi4 Frequency)')
    #plt.axvline(x=0.028)
    #plt.xlim(0,0.2)
    plt.xlim(500,2000)
    #plt.ylim(-3,2.5)
    
plt.tight_layout()
plt.legend()
plt.show()
plt.close()
In [ ]:
#Strain vs Time


color = ['C0','C1','C2','r', 'b', 'k','m','coral','forestgreen']
ls = ['-', '--', ':','-.','--',':']
j=0

fig=plt.figure(figsize=(15,6))
for i, key in enumerate(sorted(compare_nsbh.keys())):
    
    labelname = compare_nsbh[key]['label']
    
    if (compare_nsbh[key]['type']=='NSBH'):
        datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
        #filepath = os.path.join(datadir, 'Strain_l2_m2.txt')
    elif(compare_nsbh[key]['type']=='BBH'):
        datadir = compare_nsbh[key]['dir']
        
    filepath = os.path.join(datadir, 'Strain/Strain_l2_m2.txt')
    if not os.path.exists(filepath):continue
    
    massratio = int((labelname.split('-')[0]).split('Q')[1])
    if massratio==2:
        j=1
    if massratio==3:
        j=2
    if massratio==5:
        j=3
    t22, str_re22, str_im22 = np.loadtxt(filepath, unpack=True, usecols=(0,1,2))
    
    strain = np.sqrt(str_re22**2 + str_im22**2)
    max_amp_idx   = np.where(strain==np.amax(strain))
    
    
    t22_ret = t22-R
    plt.subplot(1,3,j)
    plt.plot(t22_ret - t22_ret[max_amp_idx], str_re22, c=color[i], label=labelname)
    #plt.plot(t22_ret - t22_ret[max_amp_idx], np.sqrt(str_re22**2 + str_im22**2),c=color[i])
    plt.xlabel(r'$(T - R_{ex}/M_{0}$')
    plt.ylabel(r'$r \;h(t)\;/M_{0}$')
    plt.xlim(-100, 100)
    #plt.ylim(1e-28, 1e-23)
    plt.legend()

    
plt.tight_layout()
plt.legend()
plt.grid(True)
fig.suptitle('Figure 4.1: Strain h(t) Comparison', fontsize=16, y=1.08)
plt.show()
plt.close()


#Strain vs Time


color = ['C0','C1','C2','r', 'b', 'k','m','coral','forestgreen']
ls = ['-', '--', ':','-.','--',':']


fig=plt.figure(figsize=(15,6))
for i, key in enumerate(sorted(compare_nsbh.keys())):
    
    labelname = compare_nsbh[key]['label']
    massratio = int((labelname.split('-')[0]).split('Q')[1])
    
    if (compare_nsbh[key]['type']=='NSBH'):
        datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
        #filepath = os.path.join(datadir, 'Strain_l2_m2.txt')
    elif(compare_nsbh[key]['type']=='BBH'):
        datadir = compare_nsbh[key]['dir']
        
    filepath = os.path.join(datadir, 'Strain/Strain_l2_m2.txt')
    if not os.path.exists(filepath):continue
    
    t22, str_re22, str_im22 = np.loadtxt(filepath, unpack=True, usecols=(0,1,2))
    
    strain = np.sqrt(str_re22**2 + str_im22**2)
    str_phase = np.unwrap(np.arctan2(str_im22, str_re22))
    str_freq = np.gradient(str_phase)/np.gradient(t22)/np.pi/2.
    max_amp_idx   = np.where(strain==np.amax(strain))
    M_in_s1 = (1+massratio)*1.35/s1_in_Msun
    #print(labelname,t22[np.argmax(t22>150)],str_freq[np.argmax(t22>150)], str_freq[np.argmax(t22>150)]/M_in_s1)
    
    str_physunits = (strain/0.1)*2.4*1e-22*(1+massratio)*1.35/5.
    t22_src = t22
    
    plt.subplot(121)
    #plt.plot(t22_ret - t22_ret[max_amp_idx], str_re22, c=color[i])
    plt.plot(t22 - t22[max_amp_idx], strain,c=color[i], label=labelname)
    #plt.plot((t22 - t22[max_amp_idx])*M_in_s1, str_physunits,c=color[i], label=labelname)
    plt.xlabel(r'$(T - R_{ex}/M_{0}$')
    plt.ylabel(r'$r \;h(t)\;/M_{0}$')
    plt.xlim(-100, 100)
    #plt.xlim(-0.01, 0.003)
    #plt.ylim(1e-28, 1e-23)
    plt.legend()
    plt.grid(True)

    plt.subplot(122)
    #plt.plot(t22_ret - t22_ret[max_amp_idx], str_re22, c=color[i])
    plt.plot(t22_src - t22_src[max_amp_idx],str_freq ,c=color[i], label=labelname)
    plt.xlabel(r'$(T - R_{ex})/M_{0}$')
    plt.ylabel(r'$ Freq$')

    plt.xlim(-100, 100)
    plt.ylim(0,0.1)
    plt.legend()
    plt.grid(True)

   
    
plt.tight_layout()
plt.legend()

fig.suptitle('Figure 4.2: Strain h(t) Amplitude Comparison', fontsize=16, y=1.08)
plt.show()
plt.close()
In [ ]:
#QNM Fitting - Thesis Figure 
#Psi4 plots: BBH vs NSBH
from collections import defaultdict

i=0


color = ['k','b', 'orangered','r', 'b', 'k','m','coral','forestgreen']
ls = ['-', '--', ':','-.','--',':']


lm_modes = {'2,2':(2,2), '2,1':(2,1), '3,3':(3,3), '3,2':(3,2)}#, '4,4':(4,4), '4,3':(4,3)}
y_amp_max = [1e-4, 5e-4, 1e-4,  ]
y_phase_max = [120, 150,50,250]
y_phase_min = [-50, 20,-50,50]
qnm_tau_dict = defaultdict(dict)
qnm_freq_dict = defaultdict(dict)
fig = plt.figure(figsize=(14,16))
#fig = plt.figure(figsize=(14,20))

nsbh_iter = 0
  
for lm_idx, lm_key in enumerate(sorted(lm_modes.keys())):

    lm=lm_modes[lm_key]
    #print("\nModes: l=%d, m=%d \n"%(lm[0], lm[1]))
    #print("Simulation \t QNM Tau \t QNM Freqeuncy")
    #print("-"*50)

    for i, key in enumerate(sorted(compare_nsbh.keys())):
            
        labelname = compare_nsbh[key]['label']
        if (compare_nsbh[key]['type']=='NSBH'):
            datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
            nsbh_iter+=1
        elif(compare_nsbh[key]['type']=='BBH'):
            datadir = compare_nsbh[key]['dir']
            continue
     
    
       
        
        filepath = os.path.join(datadir, 'Ylm_WEYLSCAL4::Psi4r_l%d_m%d_r70.00.asc'%(lm[0],lm[1]))
        if not os.path.exists(filepath):
            filepath = os.path.join(datadir, 'mp_Psi4r_l%d_m%d_r70.00.asc'%(lm[0],lm[1]))
        
        t22, psi4_re22, psi4_im22 = np.loadtxt(filepath, unpack=True, usecols=(0,1,2))
    
        amp = np.sqrt(psi4_re22**2 + psi4_im22**2)
        
        phase = -1*np.unwrap(np.arctan2(psi4_im22,psi4_re22))
        #Compute QNM Modes
        #For QNM - A = constant*exp(-t\Ï„) and phase = omega*t+phi_0
        
        #Find the fitting time interval
        max_amp_idx = np.argmax(amp==np.amax(amp))
        qnm_begin_idx = np.amin(np.where(t22 > t22[max_amp_idx]+40))
        qnm_end_idx = np.amin(np.where(t22 > t22[max_amp_idx]+70))
        
        
        #Create the fits 
        qnm_t_int     = t22[qnm_begin_idx:qnm_end_idx]
        qnm_amp_int   = amp[qnm_begin_idx:qnm_end_idx]
        qnm_phase_int = phase[qnm_begin_idx:qnm_end_idx]
        
        
        logamp_fit_coeff   = np.polyfit(qnm_t_int, np.log(qnm_amp_int),1)
        logamp_fit         = np.poly1d(logamp_fit_coeff)
        phase_fit_coeff    = np.polyfit(qnm_t_int, qnm_phase_int,1)
        phase_fit          = np.poly1d(phase_fit_coeff)
        logamp_fit_slope = logamp_fit_coeff[0]
        
        qnm_tau = -1./logamp_fit_slope
        qnm_freq = phase_fit_coeff[0]/2./np.pi
        
        qnm_tau_dict[lm_key][labelname] = qnm_tau
        qnm_freq_dict[lm_key][labelname] = qnm_freq
        
        #print("%10s \t %g \t %g "%(labelname, qnm_tau, qnm_freq))
        #if lm[1]==1 or lm[1]==4: continue
        #print(lm_idx, 2*lm_idx+1, 2*lm_idx+2)
        plt.subplot(4,2,2*lm_idx+1)
        amp_plot = plt.semilogy(t22-t22[max_amp_idx], amp,ls='-', c=color[i],label=labelname[:2])
        #plt.plot(t22[max_amp_idx:]-t22[max_amp_idx], amp_fit(t22[max_amp_idx:]),ls='--', c=color[i])
        #plt.plot(qnm_t_int-t22[max_amp_idx], logamp_fit(qnm_t_int),ls='-', c=color[i])
       
        
        if lm_idx==2:plt.xlabel(r'$T/M$')
        plt.ylabel(r'$|\psi_4|$')
        plt.xlim(0,130)
        #plt.ylim(1e-9,4e-4)
        plt.legend()
        plt.grid(True)
        
        
        plt.subplot(4,2,2*lm_idx+2)
        plt.plot(t22 - t22[max_amp_idx], phase,ls='-', c=color[i], label=labelname[:2])
        #plt.plot(qnm_t_int - t22[max_amp_idx], phase_fit(qnm_t_int), ls='-', c=color[i])#, label=labelname+'Fit')
        
        if lm_idx==2:plt.xlabel(r'$T/M$')
        plt.ylabel(r'$\phi$')
        
        plt.xlim(0,130)
        plt.ylim(y_phase_min[lm_idx], y_phase_max[lm_idx])
        plt.legend()
        plt.grid(True)
    
#fig.suptitle('(%s) Mode'%lm_key, fontsize=16, y=1.08)
plt.tight_layout()
#plt.savefig(os.path.join(figdir, 'QNMFitting.png'), dpi=100)

plt.show()
plt.close()
In [ ]:
# Fit QNM Omega (from Berti's datafiles) using python fitting function

def fit_omega(a_bh, f1, f2, f3):
    
    Momega_lmn = (f1 + f2*(1.-a_bh)**f3)
    return Momega_lmn

def fit_Q(a_bh, q1, q2, q3):
    
    Q_lmn = (q1 + q2*(1.-a_bh)**q3)
    return Q_lmn

lmn = '210'
j, Momega_re, Momega_img = np.loadtxt('/nethome/bkhamesra3/Desktop/rsync/l2/n1l%sm%s.dat'%(lmn[0], lmn[1]), unpack=True, usecols=(0,1,2))

tau_by_M = -1./Mtau_inv
Q = Momega_re*tau_by_M/2.0
popt_omega, pcov_omega = curve_fit(fit_omega, j, Momega_re)
popt_tau, pcov_tau = curve_fit(fit_Q, j, Q)
a = 0.6

print(1./(Q[j==a]*2/Momega_re[j==a]), 1./(fit_Q(a, *popt_tau)*2/fit_omega(a, *popt_omega)))

plt.figure(figsize=(15,6))
plt.subplot(121)
plt.plot(j, Momega_re, label='Data')
plt.plot(j, fit_omega(j, *popt_omega), label='fit - Bhavesh')
plt.plot(j, fit_omega(j, f1_dict[lmn], f2_dict[lmn],f3_dict[lmn]), label='fit -Berti Paper')
plt.xlabel('a')
plt.ylabel(r'$\omega$')
plt.legend()

plt.subplot(122)
plt.plot(j, (Momega_re - fit_omega(j, *popt_omega))/Momega_re*100, label='fit - Bhavesh')
plt.plot(j,(Momega_re -  fit_omega(j, f1_dict[lmn], f2_dict[lmn],f3_dict[lmn]))/Momega_re*100, label='fit - Berti')
plt.xlabel('a')
plt.ylabel(r'$\omega$ ')
plt.legend()
plt.show()
plt.close()


plt.figure(figsize=(15,6))
plt.subplot(121)
plt.plot(j, Q, label='Data')
plt.plot(j, fit_Q(j, *popt_tau), label='fit - Bhavesh')
plt.plot(j, fit_Q(j, q1_dict[lmn], q2_dict[lmn],q3_dict[lmn]), label='fit -Berti Paper')
plt.xlabel('a')
plt.ylabel('Q ')
plt.legend()

plt.subplot(122)
plt.plot(j, (Q - fit_Q(j, *popt_tau))/Q*100, label='fit - Bhavesh')
plt.plot(j,(Q -  fit_Q(j, q1_dict[lmn], q2_dict[lmn],q3_dict[lmn]))/Q*100, label='fit -Berti Paper')
plt.xlabel('a')
plt.ylabel('Q % error')
plt.legend()
plt.show()
plt.close()
In [ ]:
#Strain - Mode Mixing - Thesis Figure


color = ['C0','C1','C2','r', 'b', 'k','m','coral','forestgreen']
ls = ['-', '--', ':','-.','--',':']
j=0

fig=plt.figure(figsize=(15,6))
for i, key in enumerate(sorted(compare_nsbh.keys())):
    
    labelname = compare_nsbh[key]['label']
    
    if (compare_nsbh[key]['type']=='NSBH'):
        datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
        #filepath = os.path.join(datadir, 'Strain_l2_m2.txt')
    elif(compare_nsbh[key]['type']=='BBH'):
        datadir = compare_nsbh[key]['dir']
    
    if  'Q5' not in labelname: continue
    
    lm_modes = [(2,1), (2,2), (3,3), (4,4)]
    for j,lm in enumerate(lm_modes):
        filepath = os.path.join(datadir, 'Strain/Strain_l%d_m%d.txt'%(lm[0], lm[1]))
        if not os.path.exists(filepath):continue
    
        massratio = int((labelname.split('-')[0]).split('Q')[1])
       
        if 'NS' in labelname:
            plt.subplot(121)
            ls='-'
        else:
            plt.subplot(122)
            ls='--'
        t22, str_re22, str_im22 = np.loadtxt(filepath, unpack=True, usecols=(0,1,2))
    
        strain = np.sqrt(str_re22**2 + str_im22**2)
        max_amp_idx   = np.where(strain==np.amax(strain))
    
        t22_ret = t22-R
        #plt.plot(t22_ret - t22_ret[max_amp_idx], str_re22, c=color[i], label=labelname)
        plt.semilogy(t22_ret -t22_ret[max_amp_idx], np.sqrt(str_re22**2 + str_im22**2), c=color[j],ls=ls,label='(%d, %d)'%(lm[0],lm[1]))
        plt.xlabel(r'$(T - R_{ex})/M_{0}$')
        plt.ylabel(r'$r \;h(t)\;/M_{0}$')
        plt.xlim(t22_ret[0] -t22_ret[max_amp_idx]+150, 50)
        plt.ylim(1e-3, 1e0)
        plt.legend()
        plt.grid(True)
        plt.title(r'$%s - q=5$'%labelname[3:], y=1.01, fontsize=18)



plt.subplot(122)
    
lm_modes = [(2,1), (2,2), (3,3), (4,4)]
for j,lm in enumerate(lm_modes):
    filepath = os.path.join(datadir, 'Strain/Strain_l%d_m%d.txt'%(lm[0], lm[1]))
    if not os.path.exists(filepath):continue
    
    t22, str_re22, str_im22 = np.loadtxt(filepath, unpack=True, usecols=(0,1,2))

    strain = np.sqrt(str_re22**2 + str_im22**2)
    max_amp_idx   = np.where(strain==np.amax(strain))

    t22_ret = t22-R
    #plt.plot(t22_ret - t22_ret[max_amp_idx], str_re22, c=color[i], label=labelname)
    plt.semilogy(t22_ret -t22_ret[max_amp_idx], np.sqrt(str_re22**2 + str_im22**2), ls='-', c=color[j])#,label='(%d, %d)'%(lm[0],lm[1]))
    plt.xlabel(r'$(T - R_{ex})/M_{0}$')
    plt.ylabel(r'$r \;h(t)\;/M_{0}$')
    plt.xlim(t22_ret[0] -t22_ret[max_amp_idx]+150, 50)
    plt.ylim(1e-3, 1e0)
    plt.legend()
    plt.grid(True)
    #plt.title(r'$%s$'%labelname[3:], y=1.01)


plt.tight_layout()
plt.savefig(os.path.join(figdir,'Strain_ModeMixing.png'), dpi=100)
#fig.suptitle('Figure 4.1: Strain h(t) Comparison', fontsize=16, y=1.08)
plt.show()
plt.close()

Remnant BH Kicks

In [ ]:
#Remnant BH Kick
from scipy.constants import c as vc

plt.figure(figsize=(12,6))
for i, key in enumerate(sorted(compare_nsbh.keys())): #[0:3:2])):
    
    labelname = compare_nsbh[key]['label']
    if (compare_nsbh[key]['type']=='NSBH'):
        datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
        rad_arr = np.array((50,60,70,80,90,100,110,120,130,140))
    elif(compare_nsbh[key]['type']=='BBH'):
        datadir = compare_nsbh[key]['dir']
        rad_arr = np.array((50,60,70,80,90,100,110,120,130,140))
    Kmag_arr = []
    for r in rad_arr:
        filepath = os.path.join(datadir, 'psi4radAway_l6_r%d.00.asc'%r)
        if not os.path.exists(filepath):continue
        t, Erad, Kx, Ky, Kz, Edot, Pxdot, Pydot, Pzdot, Jxdot, Jydot, Jzdot, Jx, Jy, Jz\
                 = np.loadtxt(filepath, unpack=True, usecols=(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14))
        Kmag = (vc/1000)*(Kx**2 + Ky**2 + Kz**2)**0.5
        Kmag_arr.append(Kmag[-1])
        
        #print ('{:<16s} \t {:>1.2g} \t {:>1.2g} \t {:>1.2g} \t {:>1.4g}'.format(labelname,Kx[-1],Ky[-1],Kz[-1], Kmag[-1] ))

    compare_nsbh[key]['Kick'] = Kmag_arr
    compare_nsbh[key]['Det_radius'] = rad_arr
    if len(Kmag_arr)<1: continue
    plt.plot(compare_nsbh[key]['Det_radius'], compare_nsbh[key]['Kick'], label=labelname)
    
plt.xlabel('Radius/M')
plt.ylabel('Kick Velocity')
plt.legend()
plt.show()
plt.close()

Energy

In [ ]:
#Energy Plot:

i=0

plt.figure(figsize=(15,8))

color = [ 'k','b','orangered','k','b','orangered' ,'k','m','coral','forestgreen']
ls = ['-', '--', ':','-.','--',':']


for i, key in enumerate(sorted(compare_nsbh.keys())): #[0:3:2])):
    
    labelname = compare_nsbh[key]['label']
    if (compare_nsbh[key]['type']=='NSBH'):
        datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
        ls='--'
    elif(compare_nsbh[key]['type']=='BBH'):
        datadir = compare_nsbh[key]['dir']
        ls='-'
        
    filepath = os.path.join(datadir, 'psi4analysis_r70.00.asc')#'psi4radAway_l6_r70.00.asc')
    
    if not os.path.exists(filepath):
        print('File note found for {}'.format(labelname))
        
    t, Erad   = np.loadtxt(filepath, unpack=True, usecols=(0,1))

    psi4_filepath = os.path.join(datadir, 'Ylm_WEYLSCAL4::Psi4r_l2_m2_r70.00.asc')
    if not os.path.exists(psi4_filepath):
        psi4_filepath = os.path.join(datadir, 'mp_Psi4r_l2_m2_r70.00.asc')
        
    t22, psi4_re22, psi4_im22 = np.loadtxt(psi4_filepath, unpack=True, usecols=(0,1,2))
    amp = np.sqrt(psi4_re22**2 + psi4_im22**2)
    max_amp_idx = np.argmax(amp==np.amax(amp))
    
    plt.subplot(121)
    plt.plot(t-t22[max_amp_idx], Erad,ls=ls, c=color[i], label=labelname)
    plt.xlabel(r'$Time - T_{max}$')
    plt.ylabel('Radiated Energy')
    plt.legend()
    
    plt.subplot(122)
    plt.plot(t-t22[max_amp_idx], Erad,ls=ls, c=color[i], label=labelname)
    plt.xlabel(r'$Time - T_{max}$')
    plt.ylabel('Radiated Energy')
    plt.xlim(-400,200)
    plt.legend()
    
plt.tight_layout()
plt.show()
plt.close()

plt.figure(figsize=(15,8))
for i, key in enumerate(sorted(compare_nsbh.keys())): #[0:3:2])):
    
    labelname = compare_nsbh[key]['label']
    if (compare_nsbh[key]['type']=='NSBH'):
        ls='-'
        datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
    elif(compare_nsbh[key]['type']=='BBH'):
        datadir = compare_nsbh[key]['dir']
        ls='--'
        continue
    filepath = os.path.join(datadir, 'psi4radAway_l6_r70.00.asc')
    if not os.path.exists(filepath):continue
    t, Erad, Kx, Ky, Kz, Edot, Pxdot, Pydot, Pzdot, Jxdot, Jydot, Jzdot, Jx, Jy, Jz\
             = np.loadtxt(filepath, unpack=True, usecols=(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14))

    psi4_filepath = os.path.join(datadir, 'Ylm_WEYLSCAL4::Psi4r_l2_m2_r70.00.asc')
    if not os.path.exists(psi4_filepath):
        psi4_filepath = os.path.join(datadir, 'mp_Psi4r_l2_m2_r70.00.asc')
    
    t22, psi4_re22, psi4_im22 = np.loadtxt(psi4_filepath, unpack=True, usecols=(0,1,2))
    amp = np.sqrt(psi4_re22**2 + psi4_im22**2)
    max_amp_idx = np.argmax(amp==np.amax(amp))
    
    plt.subplot(121)
    plt.plot(t-t22[max_amp_idx], Edot,ls=ls, c=color[i], label=labelname)
    plt.xlabel(r'$Time - T_{max}$')
    plt.ylabel('Luminosity')
    plt.legend()
    
    plt.subplot(122)
    plt.plot(t-t22[max_amp_idx], Edot,ls=ls, c=color[i], label=labelname)
    plt.xlabel(r'$Time - T_{max}$')
    plt.ylabel('Luminosity')
    plt.xlim(-300,200)
    plt.legend()
    
plt.tight_layout()
plt.show()
plt.close()
In [ ]:
#Radiated Energy - Thesis Figure
i=0

plt.figure(figsize=(15,8))

color = ['r', 'b', 'k','m','coral','forestgreen']
ls = ['-', '--', ':','-.','--',':']


for i, key in enumerate(sorted(compare_nsbh.keys())): #[0:3:2])):
    
    labelname = compare_nsbh[key]['label']
    if (compare_nsbh[key]['type']=='NSBH'):
        datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
    elif(compare_nsbh[key]['type']=='BBH'):
        datadir = compare_nsbh[key]['dir']
        
    filepath = os.path.join(datadir, 'psi4analysis_r70.00.asc')#'psi4radAway_l6_r70.00.asc')
    
    if not os.path.exists(filepath):
        print('File note found for {}'.format(labelname))
        
    t, Erad, Kx, Ky, Kz, Edot, Pxdot, Pydot, Pzdot, Jxdot, Jydot, Jzdot, Jx, Jy, Jz\
             = np.loadtxt(filepath, unpack=True, usecols=(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14))

    psi4_filepath = os.path.join(datadir, 'Ylm_WEYLSCAL4::Psi4r_l2_m2_r70.00.asc')
    if not os.path.exists(psi4_filepath):
        psi4_filepath = os.path.join(datadir, 'mp_Psi4r_l2_m2_r70.00.asc')
        
    t22, psi4_re22, psi4_im22 = np.loadtxt(psi4_filepath, unpack=True, usecols=(0,1,2))
    amp = np.sqrt(psi4_re22**2 + psi4_im22**2)
    max_amp_idx = np.argmax(amp==np.amax(amp))
    
    
    plt.subplot(121)
    plt.plot(t-t22[max_amp_idx], Erad,ls=ls[i], c=color[i], label=labelname)
    plt.xlabel(r'$Time - T_{max}$')
    plt.ylabel('Radiated Energy')
    plt.xlim(-400,200)
    plt.legend()
    
    plt.subplot(122)
    plt.plot(t-t22[max_amp_idx], Edot,ls=ls[i], c=color[i], label=labelname)
    plt.xlabel(r'$Time - T_{max}$')
    plt.ylabel('Luminosity')
    plt.xlim(-150,100)
    plt.legend()
    
    
plt.tight_layout()
plt.savefig(os.path.join(figdir, 'RadiatedEnergy.png'), dpi=100)
plt.show()
plt.close()

Anuglar Momentum (Jz)

In [ ]:
#Angular Momentum Plots:

i=0

fig = plt.figure(figsize=(15,8))

color = ['k','b', 'orangered','m','coral','forestgreen']
ls = ['-', '--', ':','-.','--',':']

for i, key in enumerate(sorted(compare_nsbh.keys())): #[0:3:2])):
    
    labelname = compare_nsbh[key]['label']
    if (compare_nsbh[key]['type']=='NSBH'):
        datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
        ihspin_filepath = os.path.join(datadir, 'ihspin_hn_0.asc')
        ls = '-'
        c = color[i-3]
    elif(compare_nsbh[key]['type']=='BBH'):
        datadir = compare_nsbh[key]['dir']
        ls='--'
        c = color[i]
        ihspin_filepath = os.path.join(datadir, 'ihspin_hn_3.asc')
    
    filepath = os.path.join(datadir, 'psi4radAway_l6_r70.00.asc')

    t, Erad, Kx, Ky, Kz, Edot, Pxdot, Pydot, Pzdot, Jxdot, Jydot, Jzdot, Jx, Jy, Jz\
             = np.loadtxt(filepath, unpack=True, usecols=(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14))
    
    psi4_filepath = os.path.join(datadir, 'Ylm_WEYLSCAL4::Psi4r_l2_m2_r70.00.asc')
    if not os.path.exists(psi4_filepath):
        psi4_filepath = os.path.join(datadir, 'mp_Psi4r_l2_m2_r70.00.asc')
    
    t22, psi4_re22, psi4_im22 = np.loadtxt(psi4_filepath, unpack=True, usecols=(0,1,2))
    amp = np.sqrt(psi4_re22**2 + psi4_im22**2)
    max_amp_idx = np.argmax(amp==np.amax(amp))
    
    #filepath = os.path.join(datadir, 'psi4analysis_r70.00.asc')

    #t2, Erad2, Kx2, Ky2, Kz2, Edot2, Pxdot2, Pydot2, Pzdot2, Jxdot2, Jydot2, Jzdot, Jx2, Jy2, Jz2\
    #         = np.loadtxt(filepath, unpack=True, usecols=(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14))
    
    
    
    plt.subplot(121)
    plt.plot(t-t22[max_amp_idx], Jz,ls=ls, c=c, label=labelname)
    #plt.plot(t2-t22[max_amp_idx], Jz2,ls=ls[i-1], c=color[i-1], label=labelname+'psi4analysis')
    
    plt.xlabel(r'$T/M$')
    plt.ylabel(r'$J_z$')
    plt.legend()
    
    plt.subplot(122)
    plt.plot(t-t22[max_amp_idx], Jz,ls=ls, c=c, label=labelname)
    #plt.plot(t2-t22[max_amp_idx], Jz2,ls=ls[i-1], c=color[i-1], label=labelname+'psi4analysis')
    
    plt.xlabel(r'$T/M$')
    plt.ylabel(r'$J_z$')
    plt.xlim(-50,80)
    plt.legend()
  

print('\n')
#fig.suptitle('Angular Momentum', y=1.08, fontsize=20)
plt.tight_layout()
plt.savefig(os.path.join(figdir, 'AngMom_Jz_Plot1.png'))
plt.show()
plt.close()
In [ ]:
#Angular Momentum Plots:

i=0

fig = plt.figure(figsize=(15,8))

color = ['k','b', 'orangered','m','coral','forestgreen']
ls = ['-', '--', ':','-.','--',':']

for i, key in enumerate(sorted(compare_nsbh.keys())): #[0:3:2])):
    
    labelname = compare_nsbh[key]['label']
    if (compare_nsbh[key]['type']=='NSBH'):
        datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
        ihspin_filepath = os.path.join(datadir, 'ihspin_hn_0.asc')
        ls = '-'
        c = color[i-3]
    elif(compare_nsbh[key]['type']=='BBH'):
        datadir = compare_nsbh[key]['dir']
        ls='--'
        c = color[i]
        ihspin_filepath = os.path.join(datadir, 'ihspin_hn_3.asc')
    
    filepath = os.path.join(datadir, 'psi4radAway_l6_r70.00.asc')

    t, Erad, Kx, Ky, Kz, Edot, Pxdot, Pydot, Pzdot, Jxdot, Jydot, Jzdot, Jx, Jy, Jz\
             = np.loadtxt(filepath, unpack=True, usecols=(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14))
    
    t_bh, sx_bh, sy_bh, sz_bh = np.loadtxt(ihspin_filepath, unpack=True, usecols=(0,1,2,3))
    
   
    psi4_filepath = os.path.join(datadir, 'Ylm_WEYLSCAL4::Psi4r_l2_m2_r70.00.asc')
    if not os.path.exists(psi4_filepath):
        psi4_filepath = os.path.join(datadir, 'mp_Psi4r_l2_m2_r70.00.asc')
    
    t22, psi4_re22, psi4_im22 = np.loadtxt(psi4_filepath, unpack=True, usecols=(0,1,2))
    amp = np.sqrt(psi4_re22**2 + psi4_im22**2)
    max_amp_idx = np.argmax(amp==np.amax(amp))
    
    #filepath = os.path.join(datadir, 'psi4analysis_r70.00.asc')

    #t2, Erad2, Kx2, Ky2, Kz2, Edot2, Pxdot2, Pydot2, Pzdot2, Jxdot2, Jydot2, Jzdot, Jx2, Jy2, Jz2\
    #         = np.loadtxt(filepath, unpack=True, usecols=(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14))
    
    
    
    plt.subplot(121)
    plt.plot(t-t22[max_amp_idx], Jz,ls=ls, c=c, label=labelname)
    #plt.plot(t2-t22[max_amp_idx], Jz2,ls=ls[i-1], c=color[i-1], label=labelname+'psi4analysis')
    
    plt.xlabel(r'$T/M$')
    plt.ylabel(r'$J_z$')
    plt.xlim(-150,300)
    plt.legend()
  
    
    plt.subplot(122)
    if 'NSBH' in labelname:
        plt.plot(t_bh-(t22[max_amp_idx]-70), sz_bh,ls=ls, c=c, label=labelname)
    if 'BBH' in labelname:
        #plt.plot(t_bh[t_bh>(t22[max_amp_idx]-70)]-(t22[max_amp_idx]-70), sz_bh[t_bh>(t22[max_amp_idx]-70)],ls=ls, c=c, label=labelname)
        plt.axhline(y = sz_bh[-1], ls=ls, color=c, label=labelname)
    
    #plt.plot(t2-t22[max_amp_idx], Jz2,ls=ls[i-1], c=color[i-1], label=labelname+'psi4analysis')
    
    plt.xlabel(r'$(T - r)/M$')
    plt.ylabel(r'$S_z$')
    plt.xlim(-150,300)
    #plt.ylim(0.4,0.5)
    plt.legend()
    
    #print(t[np.amin(np.where(t-(t22[max_amp_idx])>=300.0))], t_bh[np.amin(np.where(t_bh-(t22[max_amp_idx])>=300.0))])
    print('%s \t \t %g \t %g \t %g'%(labelname, Jz[np.amin(np.where(t-(t22[max_amp_idx]-70)>=300.0))], \
                                sz_bh[np.amin(np.where(t_bh-(t22[max_amp_idx]-70)>=300.0))], \
                                Jz[np.amin(np.where(t-(t22[max_amp_idx])>=300.0))]+ sz_bh[np.amin(np.where(t_bh-(t22[max_amp_idx]-70)>=300.0))]))
    

print('\n')
#fig.suptitle('Angular Momentum', y=1.08, fontsize=20)
plt.tight_layout()
plt.savefig(os.path.join(figdir, 'AngMom_Jz_Plot2.png'))
plt.show()
plt.close()
In [ ]:
0.40391 - 0.406509 
In [ ]:
#Angular Momentum Plots:

i=0

color = ['k','b', 'orangered','m','coral','forestgreen']
ls = ['-', '--', ':','-.','--',':']
 
fig = plt.figure(figsize=(8,8))

for i, key in enumerate(sorted(compare_nsbh.keys())): #[0:3:2])):
    

    labelname = compare_nsbh[key]['label']
    if (compare_nsbh[key]['type']=='NSBH'):
        datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
        ls = '-'
        c = color[i]
        ihspin_filepath = os.path.join(datadir, 'ihspin_hn_0.asc')

    elif(compare_nsbh[key]['type']=='BBH'):
        continue
        datadir = compare_nsbh[key]['dir']
        ihspin_filepath = os.path.join(datadir, 'ihspin_hn_3.asc')
        ls='--'
        c = color[i]
    
    adm_filepath = os.path.join(datadir, 'adm_ejp_det_1.asc')
    R_ex_adm = 100

    t_adm, Eadm, Px_adm, Py_adm, Pz_adm, Jx_adm, Jy_adm, Jz_adm   = np.loadtxt(adm_filepath, unpack=True, usecols=(0,1,2,3,4,5,6,7))

    
    #Read Psi4, Radiated Quantitiess and Spin data
    psi4rad_filepath = os.path.join(datadir, 'psi4radAway_l6_r70.00.asc')
    R_ex_gw = np.float((os.path.basename(psi4rad_filepath).split('.asc')[0]).split('_r')[1])

    t_rad, Erad, Kx_rad, Ky_rad, Kz_rad, Edot_rad, Pxdot_rad, Pydot_rad, Pzdot_rad,\
    Jxdot_rad, Jydot_rad, Jzdot_rad, Jx_rad, Jy_rad, Jz_rad\
             = np.loadtxt(psi4rad_filepath, unpack=True, usecols=(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14))
    
    
    t_bh, sx_bh, sy_bh, sz_bh = np.loadtxt(ihspin_filepath, unpack=True, usecols=(0,1,2,3))
    
    psi4_filepath = os.path.join(datadir, 'Ylm_WEYLSCAL4::Psi4r_l2_m2_r70.00.asc')
    if not os.path.exists(psi4_filepath):
        psi4_filepath = os.path.join(datadir, 'mp_Psi4r_l2_m2_r70.00.asc')
    
    
    t22, psi4_re22, psi4_im22 = np.loadtxt(psi4_filepath, unpack=True, usecols=(0,1,2))
    amp = np.sqrt(psi4_re22**2 + psi4_im22**2)
    max_amp_idx = np.argmax(amp==np.amax(amp))
    
    #Make spline fit of each and compute errors
    Jz_adm_splinefit = splrep(t_adm[::2],Jz_adm[::2],s=0,k=3)
    Jz_rad_splinefit = splrep(t_rad[::16],Jz_rad[::16],s=0,k=3)
    sz_bh_splinefit = splrep(t_bh[::16],sz_bh[::16],s=0,k=3)
     
    Jz_adm_fit = splev(t_adm,Jz_adm_splinefit, der=0)
    Jz_rad_fit = splev(t_rad, Jz_rad_splinefit, der=0) 
    Sz_fit =  splev(t_bh, sz_bh_splinefit, der=0) 
    
    t_start_fit = 70+R_ex_adm
    Jz_adm_fiterr = (Jz_adm[t_adm>t_start_fit]-Jz_adm_fit[t_adm>t_start_fit])/Jz_adm[t_adm>t_start_fit]*100
    t_Jzadm_err = t_adm[t_adm>t_start_fit] 
    t_maxerr_Jzadm = t_Jzadm_err[Jz_adm_fiterr==np.amax(Jz_adm_fiterr)]
    t_minerr_Jzadm = t_Jzadm_err[Jz_adm_fiterr==np.amin(Jz_adm_fiterr)]
    
    
    
    t_start_fit = 70+R_ex_gw
    Jz_rad_fiterr = (Jz_rad[t_rad>t_start_fit]-Jz_rad_fit[t_rad>t_start_fit])/Jz_rad[t_rad>t_start_fit]*100
    t_Jzrad_err = t_rad[t_rad>t_start_fit] 
    t_maxerr_Jzrad = t_Jzrad_err[Jz_rad_fiterr==np.amax(Jz_rad_fiterr)]
    t_minerr_Jzrad = t_Jzrad_err[Jz_rad_fiterr==np.amin(Jz_rad_fiterr)]
    
    #t_start_fit = 70
    #Sz_fiterr = (sz_bh[t_bh>t_start_fit]-Sz_fit[t_bh>t_start_fit])/sz_bh[t_bh>t_start_fit]*100
    #t_Sz_err = t_bh[t_bh>t_start_fit]
    #t_maxerr_Sz = t_Sz_err[Sz_fiterr==np.amax(Sz_fiterr)]
    #t_minerr_Sz = t_Sz_err[Sz_fiterr==np.amin(Sz_fiterr)]
    
    print(labelname)
    print("ADM Ang Mom Fitting Error: Max = %g %% at t = %g, Min = %g%% at t=%g"%(np.amax(Jz_adm_fiterr), t_maxerr_Jzadm,\
                                                                         np.amin(Jz_adm_fiterr), t_minerr_Jzadm)) 
    
    print("Rad Ang Mom Fitting Error: Max = %g %% at t = %g, Min = %g%% at t=%g"%(np.amax(Jz_fiterr), t_maxerr_Jz,\
                                                                         np.amin(Jz_fiterr), t_minerr_Jz)) 
    print("Spin Fitting Error:Max = %g %% at t = %g, Min = %g %% at t=%g"%(np.amax(Sz_fiterr), t_maxerr_Sz,
                                                                          np.amin(Sz_fiterr), t_minerr_Sz))
    
    

    #Make spline fit of each in COM frame (shift time to COM frame) 
    Jz_adm_splinefit = splrep(t_adm[::2]-R_ex_adm,Jz_adm[::2],s=0,k=3)
    Jz_rad_splinefit = splrep(t_rad[::16]-R_ex_gw,Jz_rad[::16],s=0,k=3)
    sz_bh_splinefit = splrep(t_bh[::16],sz_bh[::16],s=0,k=3)
    
    #Compute fits for common time
    Jz_adm_fit = splev(t_bh,Jz_adm_splinefit, der=0)
    Jz_rad_fit = splev(t_bh, Jz_rad_splinefit, der=0) 
    Sz_fit =  splev(t_bh, sz_bh_splinefit, der=0) 
    
    Jz_matter = Jz_adm_fit - (Jz_rad_fit + Sz_fit)
    #t_common = np.intersect1d(np.around((t_bh + r_ex), decimals=1), np.around(t, decimals=1))
    #idx_rad_gw = np.isin(np.around(t, decimals=1), t_common)
    #idx_bh = np.isin(np.around(t_bh+r_ex, decimals=1), t_common)
    #
    #print(len(t_common), len(t_bh[idx_bh]) , len(t[idx_rad_gw]))
   
    #if 'Q5' not in labelname: continue
    #plt.subplot(121)
    plt.plot(t_bh - (t22[max_amp_idx] -70), Jz_matter,ls='-', c=c, label=labelname)
    #plt.plot(t_adm , Jz_adm_fit,ls='--', c=color[i-1], label=labelname+'Fitting')
    
    plt.xlabel(r'$T /M $')
    plt.ylabel(r'$J^{remining}_z$')
    plt.ylim(0,1)
    plt.xlim(-500,500)
    plt.legend()
    
    #plt.subplot(122)
    #plt.plot(t_Jzadm_err , Jz_adm_fiterr,ls='-', c=c, label=labelname)
    #
    #plt.xlabel(r'$Time - T_{max}$')
    #plt.ylabel(r'$\Delta J_ADM(\%)$')
    ##plt.xlim(-50,100)
    #plt.legend()
  
    
    
#fig.suptitle('Angular Momentum', y=1.08, fontsize=20)
plt.tight_layout()
plt.savefig(os.path.join(figdir, 'AngMom_Jz_Matter_Plot3.png'))
plt.show()
plt.close()
    
    
In [ ]:
#BH Irreducible Mass Plots - Comparison of mass  and spin of final black hole at 300M after merger
import matplotlib.gridspec as gridspec

color = ['k','b', 'orangered','m','coral','forestgreen']
ls = ['-', '--', ':','-.','--',':']

bbh_irrmass_dict = {}
bbh_hznmass_dict = {}

print("System \t \t Init. Irr Mass    Final Irr Mass \t NSBH - BBH \t Final Spin(a) \tFinal Hzn Mass  NSBH - BBH \n ")
print('-'*115)


fig = plt.figure(constrained_layout=True, figsize=(16,8))
gs = gridspec.GridSpec(ncols=2, nrows=2, figure=fig)

ax1 = fig.add_subplot(gs[0, 0])
ax2 = fig.add_subplot(gs[0, 1])
ax3 = fig.add_subplot(gs[1,:])
for i, key in enumerate(sorted(compare_nsbh.keys())):
    
    
    labelname = compare_nsbh[key]['label']
    massratio = str((labelname.split('Q')[1]).split('-')[0])
    if (compare_nsbh[key]['type']=='BBH'):
    
        datadir = compare_nsbh[key]['dir']
        t_merge = compare_nsbh[key]['t_merger_src']
        t_comp = t_merge + 500.
        
        bhdiag1 = os.path.join(datadir,'BH_diagnostics.ah1.gp')
        t_ah1, area_ah1,irrmass_ah1, arealrad_ah1 = np.loadtxt(bhdiag1, unpack=True, usecols=(1,25,26,27))
        
        ihspin1 = os.path.join(datadir,'ihspin_hn_0.asc')
        t_ih1, sx1, sy1, sz1 = np.loadtxt(ihspin1, unpack=True, usecols=(0,1,2,3))
        
        bbh_irrmass_bh1 = irrmass_ah1[0]
        bbh_spinmag_bh1 = (sx1**2 + sy1**2  + sz1**2)**0.5
        bbh_hznmass_bh1 = (bbh_irrmass_bh1**2 + 0.25*bbh_spinmag_bh1[0]**2/bbh_irrmass_bh1**2)**0.5
        
        
        bhdiag3 = os.path.join(datadir,'BH_diagnostics.ah3.gp')
        t_ah3, area_ah3,irrmass_ah3, arealrad_ah3 = np.loadtxt(bhdiag3, unpack=True, usecols=(1,25,26,27))
        
        ihspin3 = os.path.join(datadir,'ihspin_hn_2.asc')
        t_ih3, sx3, sy3, sz3 = np.loadtxt(ihspin3, unpack=True, usecols=(0,1,2,3))
        
        #plt.plot(t_ah1-t_merge, irrmass_ah1-irrmass_ah1[0], c=color[i], ls=ls[i],label=labelname)
        #plt.plot(t_ah3-t_merge, irrmass_ah3-irrmass_ah1[0], c=color[i], ls=ls[i],label=labelname)
        
        #if 'Q2' in labelname:
        #    ax1.plot(t_ah1[:-100]-t_merge, irrmass_ah1[:-100], c=color[i], ls=ls[i],label=labelname)
        #    ax1.plot(t_ah3[530:]-t_merge, irrmass_ah3[530:], c=color[i], ls=ls[i])#,label=labelname)
        #else:
        #    ax1.plot(t_ah1-t_merge, irrmass_ah1, c=color[i], ls=ls[i],label=labelname)
        #    ax1.plot(t_ah3-t_merge, irrmass_ah3, c=color[i], ls=ls[i])#,label=labelname)
        ax1.axhline(y=irrmass_ah3[-1],c=color[i], ls='--')
        #ax2.plot(t_ah1-t_merge, np.gradient(irrmass_ah1)/np.gradient(t_ah1), c=color[i], ls=ls[i],label=labelname)
        #ax2.plot(t_ah3-t_merge, np.gradient(irrmass_ah3)/np.gradient(t_ah3), c=color[i], ls=ls[i])#,label=labelname)
    
    
        #Set all the simulations to same final time before comparing mass and spins
        area_ah3       = area_ah3    [t_ah3 < t_comp]
        irrmass_ah3    = irrmass_ah3 [t_ah3 < t_comp]
        arealrad_ah3   = arealrad_ah3[t_ah3 < t_comp]
        t_ah3          = t_ah3       [t_ah3 < t_comp]
        
        sx3   = sx3  [t_ih3 < t_comp]
        sy3   = sy3  [t_ih3 < t_comp]
        sz3   = sz3  [t_ih3 < t_comp]
        t_ih3 = t_ih3[t_ih3 < t_comp]
        
        
        bbh_irrmass_bh3 = irrmass_ah3[-1]
        bbh_spinmag_bh3 = (sx3**2 + sy3**2  + sz3**2)**0.5
        bbh_hznmass_bh3 = (bbh_irrmass_bh3**2 + 0.25*bbh_spinmag_bh3[-1]**2/bbh_irrmass_bh3**2)**0.5
        bbh_dimspin_bh3 = bbh_spinmag_bh3[-1]/bbh_hznmass_bh3**2
        
        bbh_irrmass_dict[massratio]  = bbh_irrmass_bh3
        bbh_hznmass_dict[massratio]  = bbh_hznmass_bh3
        print("%-10s \t %g \t  %g \t\t %g \t\t %g \t %g \t %g "%(labelname,bbh_irrmass_bh1, bbh_irrmass_bh3, \
                                            bbh_irrmass_bh3 - bbh_irrmass_bh3, bbh_dimspin_bh3,\
                                            bbh_hznmass_bh3,  bbh_hznmass_bh3 - bbh_hznmass_bh3))
        
        
    elif(compare_nsbh[key]['type']=='NSBH'):
        
        datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
        t_merge = compare_nsbh[key]['t_merger_src']
        t_comp = t_merge + 500.
        
        
        bhdiag = os.path.join(datadir,'BH_diagnostics.ah1.gp')
        t_ah, area_ah,irrmass_ah, arealrad_ah = np.loadtxt(bhdiag, unpack=True, usecols=(1,25,26,27))
        
        ihspin = os.path.join(datadir,'ihspin_hn_0.asc')
        t_ih, sx, sy, sz = np.loadtxt(ihspin, unpack=True, usecols=(0,1,2,3))
     
        if 'Q5' in labelname:irrmass_splinefit = splrep(t_ah[::16],irrmass_ah[::16],s=0)
        else :irrmass_splinefit = splrep(t_ah[::128],irrmass_ah[::128],s=0)
    
        #Compute mass growth rate on every 8th iteration to reduce the high frequency noise in the derivatives 
        #Truncation and round off errors leads to high frequency noise which gets amplified on computing the derivatives
        #Computing the derivative on a larger time step averages out these effect to zero)
       
        mass_growth_rate_fit =  splev(t_ah, irrmass_splinefit, der=1) 
        mass_growth_rate = np.gradient(irrmass_ah[::32])/np.gradient(t_ah[::32])
        #plt.plot(t_ah-t_merge, irrmass_ah-irrmass_ah[0], c=color[i], ls=ls[i],label=labelname)
        ax1.plot(t_ah-t_merge, irrmass_ah, c=color[i-3], ls='-',label=labelname[:2])
        #ax1.plot(t_ah-t_merge, splev(t_ah, irrmass_splinefit), color[i], ls='--')
        ax2.semilogy((t_ah-t_merge), mass_growth_rate_fit, c=color[i-3], ls='-',label=labelname[:2])
        #ax2.semilogy((t_ah-t_merge)[::32], mass_growth_rate, c=color[i], ls='--',label=labelname[:2])
        
        
        
        #Set all the simulations to same final time before comparing mass and spins
        area_ah       = area_ah    [t_ah<t_comp]
        irrmass_ah    = irrmass_ah [t_ah<t_comp]
        arealrad_ah   = arealrad_ah[t_ah<t_comp]
        t_ah          = t_ah       [t_ah<t_comp]
       
        sx   = sx  [t_ih<t_comp]
        sy   = sy  [t_ih<t_comp]
        sz   = sz  [t_ih<t_comp]
        t_ih = t_ih[t_ih<t_comp]
        
        
        
        nsbh_irrmass_bh1 = irrmass_ah[0]
        nsbh_spinmag_bh1 = (sx**2 + sy**2  + sz**2)**0.5
        nsbh_hznmass_bh1 = (nsbh_irrmass_bh1**2 + 0.25*nsbh_spinmag_bh1[0]**2/nsbh_irrmass_bh1**2)**0.5
        
        nsbh_irrmass_bh3 = irrmass_ah[-1]
        nsbh_spinmag_bh3 = (sx**2 + sy**2  + sz**2)**0.5
        nsbh_hznmass_bh3 = (nsbh_irrmass_bh3**2 + 0.25*nsbh_spinmag_bh3[-1]**2/nsbh_irrmass_bh3**2)**0.5
        nsbh_dimspin_bh3 = nsbh_spinmag_bh3[-1]/nsbh_hznmass_bh3**2
            
        print("%-1s \t %g \t  %g \t\t %g \t %g \t %g \t %g"%(labelname, nsbh_irrmass_bh1,nsbh_irrmass_bh3, \
                                                  (nsbh_irrmass_bh3-bbh_irrmass_dict[massratio])/bbh_irrmass_dict[massratio], nsbh_dimspin_bh3,\
                                                  nsbh_hznmass_bh3,   (nsbh_hznmass_bh3-bbh_hznmass_dict[massratio])/bbh_hznmass_dict[massratio]))
       
        
ax1.set_xlabel(r'$T/M$', fontsize=20)
ax1.set_ylabel(r'$ M_h$', fontsize=20)
ax2.set_xlabel(r'$T/M$', fontsize=20)
ax2.set_ylabel(r'$\dot{M}_h$', fontsize=20)

ax1.set_xlim(-200,300)
ax2.set_xlim(-50,100)
#ax1.set_ylim(-2,5)
ax2.set_ylim(2e-5,0.1)

ax1.legend(loc=4)
ax2.legend()
#plt.axhline(y=irrmass_ah[0], ls='--', c='k')
#plt.ylim(0.65, 0.92)

ax1.grid(True)
ax2.grid(True)
plt.tight_layout()
#plt.savefig(os.path.join(figdir, 'BH_MassPlots.png'), dpi=100)
plt.show()
plt.close()
In [ ]:
#Rate of change of Angular Momentum 
plt.figure(figsize=(15,8))
for i, key in enumerate(sorted(compare_nsbh.keys())): #[0:3:2])):
    
    labelname = compare_nsbh[key]['label']
    if (compare_nsbh[key]['type']=='NSBH'):
        datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
        ls = '-'
        c = color[i-3]
    elif(compare_nsbh[key]['type']=='BBH'):
        datadir = compare_nsbh[key]['dir']
        ls = '--'
        c = color[i]
    
    filepath = os.path.join(datadir, 'psi4radAway_l6_r70.00.asc')

    t, Erad, Kx, Ky, Kz, Edot, Pxdot, Pydot, Pzdot, Jxdot, Jydot, Jzdot, Jx, Jy, Jz\
             = np.loadtxt(filepath, unpack=True, usecols=(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14))

    psi4_filepath = os.path.join(datadir, 'Ylm_WEYLSCAL4::Psi4r_l2_m2_r70.00.asc')
    if not os.path.exists(psi4_filepath):
        psi4_filepath = os.path.join(datadir, 'mp_Psi4r_l2_m2_r70.00.asc')
    
    t22, psi4_re22, psi4_im22 = np.loadtxt(psi4_filepath, unpack=True, usecols=(0,1,2))
    amp = np.sqrt(psi4_re22**2 + psi4_im22**2)
    max_amp_idx = np.argmax(amp==np.amax(amp))
    

    plt.subplot(121)
    plt.plot(t-t22[max_amp_idx], Jzdot,ls=ls, c=c, label=labelname)
    plt.xlabel(r'$Time - T_{max}$')
    plt.ylabel(r'$dJ_z/dt$')
    plt.legend()
    
    plt.subplot(122)
    plt.plot(t-t22[max_amp_idx], Jzdot,ls=ls, c=c, label=labelname)
    plt.xlabel(r'$Time - T_{max}$')
    plt.ylabel(r'$dJ_z/dt$')
    plt.xlim(-50,80)
    plt.legend()
    
plt.tight_layout()
plt.show()
plt.close()
In [ ]:
#Remnant Kicks:

i=0

plt.figure(figsize=(15,8))

color = ['r', 'b', 'k','m','coral','forestgreen']
ls = ['-', '--', ':','-.','--',':']


for i, key in enumerate(sorted(compare_nsbh.keys())): #[0:3:2])):
    
    labelname = compare_nsbh[key]['label']
    if (compare_nsbh[key]['type']=='NSBH'):
        datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
    elif(compare_nsbh[key]['type']=='BBH'):
        datadir = compare_nsbh[key]['dir']
    
    filepath = os.path.join(datadir, 'psi4radAway_l6_r70.00.asc')

    t, Erad, Kx, Ky, Kz, Edot, Pxdot, Pydot, Pzdot, Jxdot, Jydot, Jzdot, Jx, Jy, Jz\
             = np.loadtxt(filepath, unpack=True, usecols=(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14))

    psi4_filepath = os.path.join(datadir, 'Ylm_WEYLSCAL4::Psi4r_l2_m2_r70.00.asc')
    if not os.path.exists(psi4_filepath):
        psi4_filepath = os.path.join(datadir, 'mp_Psi4r_l2_m2_r70.00.asc')
    
    t22, psi4_re22, psi4_im22 = np.loadtxt(psi4_filepath, unpack=True, usecols=(0,1,2))
    amp = np.sqrt(psi4_re22**2 + psi4_im22**2)
    max_amp_idx = np.argmax(amp==np.amax(amp))
    
    plt.subplot(121)
    plt.plot(t-t22[max_amp_idx], Kx,ls=ls[i], c=color[i], label=labelname)
    
    plt.xlabel(r'$Time - T_{max}$')
    plt.ylabel('Kx')
    plt.xlim(-50,80)
    plt.legend()
    
    plt.subplot(122)
    plt.plot(t-t22[max_amp_idx], Ky,ls=ls[i], c=color[i], label=labelname)
    plt.xlabel(r'$Time - T_{max}$')
    plt.ylabel('Ky')
    plt.xlim(-50,80)
    plt.legend()
    
plt.tight_layout()
plt.show()
plt.close()

plt.figure(figsize=(15,8))
for i, key in enumerate(sorted(compare_nsbh.keys())): #[0:3:2])):
    
    labelname = compare_nsbh[key]['label']
    if (compare_nsbh[key]['type']=='NSBH'):
        datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
    elif(compare_nsbh[key]['type']=='BBH'):
        datadir = compare_nsbh[key]['dir']
    
    filepath = os.path.join(datadir, 'psi4radAway_l6_r70.00.asc')

    t, Erad, Kx, Ky, Kz, Edot, Pxdot, Pydot, Pzdot, Jxdot, Jydot, Jzdot, Jx, Jy, Jz\
             = np.loadtxt(filepath, unpack=True, usecols=(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14))

    psi4_filepath = os.path.join(datadir, 'Ylm_WEYLSCAL4::Psi4r_l2_m2_r70.00.asc')
    if not os.path.exists(psi4_filepath):
        psi4_filepath = os.path.join(datadir, 'mp_Psi4r_l2_m2_r70.00.asc')
    
    t22, psi4_re22, psi4_im22 = np.loadtxt(psi4_filepath, unpack=True, usecols=(0,1,2))
    amp = np.sqrt(psi4_re22**2 + psi4_im22**2)
    max_amp_idx = np.argmax(amp==np.amax(amp))
    
    plt.subplot(121)
    plt.plot(t-t22[max_amp_idx], Pxdot,ls=ls[i], c=color[i], label=labelname)
    plt.xlabel(r'$Time - T_{max}$')
    plt.ylabel(r'$dP_x/dt$')
    plt.xlim(-50,80)
    plt.legend()
    
    
    plt.subplot(122)
    plt.plot(t-t22[max_amp_idx], Pydot,ls=ls[i], c=color[i], label=labelname)
    plt.xlabel(r'$Time - T_{max}$')
    plt.ylabel(r'$dP_y/dt$')
    plt.xlim(-50,80)
    plt.legend()
    
plt.tight_layout()
plt.show()
plt.close()

ADM Quantities

In [ ]:
#ADM Energy Plot:

i=0

plt.figure(figsize=(15,8))

color = ['r', 'b', 'k','m','coral','forestgreen']
ls = ['-', '--', ':','-.','--',':']


for i, key in enumerate(sorted(compare_nsbh.keys())):
    
    labelname = compare_nsbh[key]['label']
    if (compare_nsbh[key]['type']=='NSBH'):
        datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
    elif(compare_nsbh[key]['type']=='BBH'):
        continue
        datadir = compare_nsbh[key]['dir']
    
    filepath = os.path.join(datadir, 'adm_ejp_det_3.asc')

    t, Eadm, Px_adm, Py_adm, Pz_adm, Jx_adm, Jy_adm, Jz_adm   = np.loadtxt(filepath, unpack=True, usecols=(0,1,2,3,4,5,6,7))

    psi4_filepath = os.path.join(datadir, 'Ylm_WEYLSCAL4::Psi4r_l2_m2_r70.00.asc')
    if not os.path.exists(psi4_filepath):
        psi4_filepath = os.path.join(datadir, 'mp_Psi4r_l2_m2_r70.00.asc')
    
    t22, psi4_re22, psi4_im22 = np.loadtxt(psi4_filepath, unpack=True, usecols=(0,1,2))
    amp = np.sqrt(psi4_re22**2 + psi4_im22**2)
    max_amp_idx = np.argmax(amp==np.amax(amp))
    
    plt.subplot(121)
    plt.plot(t-t22[max_amp_idx], Eadm,ls=ls[i], c=color[i], label=labelname)
    plt.xlabel(r'$Time - T_{max}$')
    plt.ylabel('ADM Energy')
    plt.legend()
    
    plt.subplot(122)
    plt.plot(t-t22[max_amp_idx], Eadm,ls=ls[i], c=color[i], label=labelname)
    plt.xlabel(r'$Time - T_{max}$')
    plt.ylabel('ADM Energy')
    plt.xlim(-50,80)
    plt.legend()
    
plt.tight_layout()
plt.show()
plt.close()

  
In [ ]:
#ADM Angular momentum Plot:

i=0

plt.figure(figsize=(10,8))

color = ['r', 'b', 'k','m','coral','forestgreen']
ls = ['-', '--', ':','-.','--',':']


for i, key in enumerate(sorted(compare_nsbh.keys())):
    
    labelname = compare_nsbh[key]['label']
    if (compare_nsbh[key]['type']=='NSBH'):
        datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
    elif(compare_nsbh[key]['type']=='BBH'):
        datadir = compare_nsbh[key]['dir']
        continue
    
    filepath = os.path.join(datadir, 'adm_ejp_det_2.asc')

    t, Eadm, Px_adm, Py_adm, Pz_adm, Jx_adm, Jy_adm, Jz_adm   = np.loadtxt(filepath, unpack=True, usecols=(0,1,2,3,4,5,6,7))

    psi4_filepath = os.path.join(datadir, 'Ylm_WEYLSCAL4::Psi4r_l2_m2_r70.00.asc')
    if not os.path.exists(psi4_filepath):
        psi4_filepath = os.path.join(datadir, 'mp_Psi4r_l2_m2_r70.00.asc')
    
    t22, psi4_re22, psi4_im22 = np.loadtxt(psi4_filepath, unpack=True, usecols=(0,1,2))
    amp = np.sqrt(psi4_re22**2 + psi4_im22**2)
    max_amp_idx = np.argmax(amp==np.amax(amp))
    
    #plt.subplot(231)
    #plt.plot(t-t22[max_amp_idx], Jx_adm,ls=ls[i], c=color[i], label=labelname)
    #plt.xlabel(r'$Time - T_{max}$')
    #plt.ylabel('$J_x^{ADM}$')
    #plt.legend()
    #
    #plt.subplot(232)
    #plt.plot(t-t22[max_amp_idx], Jx_adm,ls=ls[i], c=color[i], label=labelname)
    #plt.xlabel(r'$Time - T_{max}$')
    #plt.ylabel('$J_x^{ADM}$')
    #plt.xlim(-50,80)
    #plt.legend()
    #
    #plt.subplot(233)
    #plt.plot(t-t22[max_amp_idx], Jy_adm,ls=ls[i], c=color[i], label=labelname)
    #plt.xlabel(r'$Time - T_{max}$')
    #plt.ylabel('$J_y^{ADM}$')
    #plt.legend()
    #
    #plt.subplot(234)
    #plt.plot(t-t22[max_amp_idx], Jy_adm,ls=ls[i], c=color[i], label=labelname)
    #plt.xlabel(r'$Time - T_{max}$')
    #plt.ylabel('$J_y^{ADM}$')
    #plt.xlim(-50,80)
    #plt.legend()
    
    #plt.subplot(121)
    plt.plot(t-t22[max_amp_idx], Jz_adm,ls=ls[i], c=color[i], label=labelname)
    plt.xlabel(r'$T/M$')
    plt.ylabel('$J_z^{ADM}$')
    plt.legend()
    
    #plt.subplot(122)
    #plt.plot(t-t22[max_amp_idx], Jz_adm,ls=ls[i], c=color[i], label=labelname)
    #plt.xlabel(r'$Time - T_{max}$')
    #plt.ylabel('$J_z^{ADM}$')
    #plt.xlim(-500,500)
    #plt.legend()
    
    
plt.tight_layout()
plt.savefig(os.path.join(figdir, 'ADM_AngMom_Jz.png'))
plt.show()
plt.close()

Density Oscillations

The star is initially not quite stable with strong oscillations in cenral density. The density changes by >25% compared to the initial value. The oscillations do not reduce much continuing all the way to the merger. However, from the above analysis, it appears that these oscillations do not appear to affect the dynamics and gravitational wave signature significantly to the first order. A more rigorous comparison study to test the effects of stellar oscillations is required.

IGR vs Whisky - In both the simulations, the oscillations in the stellar density are of the same order and frequency. However, near the merger, the disruption of star and process of formation of accretion disks have significant variations. But, once the disk starts to settle, the density profile looks quite similar (check the density snapshots) in both the codes with final densities ~$10^{13}-10^{14} g/cm^3$.

The complete density and velocity evolutions can be found here

SH_P1 - Density evolutions not provided in corrected paper for q=2

Some Tests to try -

  • Compare Boosted TOV with Moving TOV obtained using BowenID
  • Do a full numerical comparison with actual simulation dataset from another independent code such as LORENE.
In [ ]:
#Density Plots

i=0
filename = 'MinSearch0.asc'
fig, ax1 = plt.subplots(1,1,figsize=(15,8))

#left, bottom, width, height = [0.62, 0.57, 0.35, 0.4]    
#ax2 = fig.add_axes([left, bottom, width, height])
color = ['r', 'b', 'k','m','coral','forestgreen']
ls = ['-', '--', ':','-.','--',':']

print("Compactness \t Oscillation Amplitude \t Relative change (%)")
for i, direc in enumerate(sorted(compare_nsbh.keys())):
    
    if compare_nsbh[direc]['type']=='BBH':
        continue
    labelname = compare_nsbh[direc]['label']
    datadir = os.path.join(compare_nsbh[direc]['dir'], 'Summary/Data')
    filepath = os.path.join(datadir, filename)
    t, max_rho = np.loadtxt(filepath, unpack=True, usecols=(1,5))
    
    idx_stablestar = np.intersect1d(np.where(t>100), np.where(t<120))
    err = (np.amax(max_rho[idx_stablestar])-np.amin(max_rho[idx_stablestar]))/np.amin(max_rho[idx_stablestar])
    if 'HR' not in labelname:
        print("%s \t \t  %g \t \t %g"%(labelname, np.amax(err*np.amin(max_rho[idx_stablestar])), np.amax(err)*100))
    
    
    ax1.plot(t, max_rho/max_rho[0], ls=ls[i],c=color[i], label=labelname)
    ax1.set_xlabel('Time')
    ax1.set_ylabel(r'$\rho_{max}/\rho_{init}$')
    ax1.legend()
    ax1.set_ylim(0.8, 1.5)
    ax1.set_xlim(200,400)
    ax1.axvline(x=230)
    ax1.axvline(x=270)
    
   
    #ax2.plot(t, max_rho/max_rho[0], ls=ls[i],c=c[i])
    #ax2.set_xlim(0, 100)
    #ax2.set_ylim(0.9, 1.3)
    #ax2.xaxis.tick_top()
    #plt.xticks([])
    #plt.yticks([])
    
    

    i+=1  
plt.tight_layout()    
plt.ylim(0,2)
#plt.xlim(0,100)
#plt.savefig()
plt.show()
plt.close()
  
    
In [ ]:
#Combine density and velocity h5 files

import subprocess
for i, key in enumerate(sorted(compare_nsbh.keys())):
    
    labelname = compare_nsbh[key]['label']
    if (compare_nsbh[key]['type']=='NSBH'):
        datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
    else:
        continue
        
    vel_h5 = sorted(glob.glob(os.path.join(datadir, 'vel???.xy.h5')))
    if len(vel_h5)<1:
        velx_h5 = [os.path.join(datadir, 'velx.xy.h5')]
        vely_h5 = [os.path.join(datadir, 'vely.xy.h5')]
        velz_h5 = [os.path.join(datadir, 'velz.xy.h5')]
        vel_h5 = velx_h5 + vely_h5 + velz_h5
        

    hdf5_merge = ["/nethome/bkhamesra3/Cactus_branches/Cactus_2019/exe/ETK/hdf5_merge"]
    rho_h5 = [os.path.join(datadir, 'rho.xy.h5')]
    h5_files = rho_h5+vel_h5
    process_string = hdf5_merge+h5_files+[os.path.join(datadir,'rho_vel.xy.h5')]

    #print process_string
    combine_h5 = subprocess.check_output(process_string) #, stdout=subprocess.PIPE)
    
In [ ]:
#Density snapshots with velocity profile

for i, key in enumerate(sorted(compare_nsbh.keys())):
    
    labelname = compare_nsbh[key]['label']
    if (compare_nsbh[key]['type']=='NSBH'):
        datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
    else:
        continue
    
    print("%s: Density - Veocity Profile Snapshots"%labelname)
     
    mpl.rc('axes', labelsize=16, grid=False)
    varname = 'density'
    filepath = os.path.join(datadir, 'rho_vel.xy.h5')
    zoom=1
    varlim=None
    annotate_grids=False                
    annotate_sphere=True
    
    bhdiag = os.path.join(datadir,'BH_diagnostics.ah1.gp')
    if os.path.exists(bhdiag):
        it_bh,time_bh,x_bh,y_bh,z_bh, rad_bh = np.loadtxt(bhdiag, unpack=True, usecols=(0,1,2,3,4,7))
        

    ds_allit = yt.load(os.path.join(datadir, 'rho_vel.xy.h5'))
    iterations = ds_allit.get_all_iterations()
    
    #Only look at merger snapshots
    #iterations = iterations[int(3*len(iterations)/4):]
    
    freq = int(len(iterations)/3)
  
    if labelname=='NSBH-Whisky': init_idx = 2
    else: init_idx=1
        
    for j in iterations:
        
        if os.path.exists(bhdiag): 
            time_j = time_bh[it_bh==j]
            
            xj,yj,zj = x_bh[it_bh==j],y_bh[it_bh==j], z_bh[it_bh==j]
            rad_j = rad_bh[it_bh==j]
        
        ds = yt.load(filepath, iteration=j)
        
        min_var, max_var = ds.all_data().quantities.extrema(varname)
        
        var = yt.SlicePlot(ds, 'z', varname,  axes_unit='code_length', origin='native', width=60, center = [0,0,0]) #center=[xnsj, ynsj, znsj])
        
        if annotate_grids:  var.annotate_grids()
        
        
        var.set_cmap(varname, 'dusk') #'viridis')
        var.annotate_velocity(plot_args={"headwidth": 2,'color':'white'})
        
        if len(time_j)>0:var.annotate_text((0.45,0.9), 't = %.2f'%time_j, coord_system='figure', text_args={'color':'white'})
        else:var.annotate_text((0.45,0.9), 'it = %d'%j, coord_system='figure', text_args={'color':'white'})
        
        if os.path.exists(bhdiag) and annotate_sphere:
            try:
                    var.annotate_sphere([xj,yj,zj], radius=(rad_j, 'code_length'), circle_args={'facecolor':'black','fill':0, 'edgecolor':'white'})
            except YTPlotCallbackError:
                print("Annotation Failed")
       
 
        try:
            var.zoom((zoom))
            var.set_font_size(20)
            var.set_xlabel('x (M)')
            var.set_ylabel('y (M)')
            if varname=='density':
                var.set_zlim('density', 1.0e9, 1e16) #max_var)
            elif varlim!=None:
                var.set_zlim(varname, varlim[0], varlim[1])
            else:
                var.set_zlim(varname, 5e-2, max_var)
        
            # This forces the ProjectionPlot to redraw itself on the AxesGrid axes.
            plot = var.plots[varname]
        
            cbar = plot.cb
        
            # Finally, this actually redraws the plot.
            var._setup_plots()
        
            if varname=='density':
            
                tick_labels = []
                ticks=10**(np.arange(12,19))
                for i in range(12,19):
                    tick_labels.append(r"$10^{%d}$"%(i))
                
           # cbar.set_ticks(ticks.tolist())    
           # cbar.set_ticklabels(tick_labels)  
        
        # Finally, this actually redraws the plot.
        #var._setup_plots()
        except  :
            print("Annotation Failed")
        # Ensure the colorbar limits match for all plots
        
        var.show()


    
In [ ]:
#BH Horizon Mass Plots - Comparison of mass  and spin of final black hole at 300M after merger
color = ['r', 'b', 'k','m','coral','forestgreen']
ls = ['-', '--', ':','-.','--',':']

bbh_irrmass_dict = {}
bbh_hznmass_dict = {}

print("System \t \t Init. Irr Mass    Final Irr Mass \t NSBH - BBH \t Final Spin(a) \tFinal Hzn Mass  NSBH - BBH \n ")
print('-'*115)


fig, (ax1, ax2) = plt.subplots(1,2,figsize=(15,6))

for i, key in enumerate(sorted(compare_nsbh.keys())):
    
    
    labelname = compare_nsbh[key]['label']
    massratio = str((labelname.split('Q')[1]).split('-')[0])
    if (compare_nsbh[key]['type']=='BBH'):
        datadir = compare_nsbh[key]['dir']
        t_merge = compare_nsbh[key]['t_merger_src']
        t_comp = t_merge + 300.
        
        bhdiag1 = os.path.join(datadir,'BH_diagnostics.ah1.gp')
        t_ah1, area_ah1,irrmass_ah1, arealrad_ah1 = np.loadtxt(bhdiag1, unpack=True, usecols=(1,25,26,27))
        
        ihspin1 = os.path.join(datadir,'ihspin_hn_0.asc')
        t_ih1, sx1, sy1, sz1 = np.loadtxt(ihspin1, unpack=True, usecols=(0,1,2,3))
        
        bbh_irrmass_bh1 = irrmass_ah1[0]
        bbh_spinmag_bh1 = (sx1**2 + sy1**2  + sz1**2)**0.5
        bbh_hznmass_bh1 = (bbh_irrmass_bh1**2 + 0.25*bbh_spinmag_bh1[0]**2/bbh_irrmass_bh1**2)**0.5
        
        
        bhdiag3 = os.path.join(datadir,'BH_diagnostics.ah3.gp')
        t_ah3, area_ah3,irrmass_ah3, arealrad_ah3 = np.loadtxt(bhdiag3, unpack=True, usecols=(1,25,26,27))
        
        ihspin3 = os.path.join(datadir,'ihspin_hn_2.asc')
        t_ih3, sx3, sy3, sz3 = np.loadtxt(ihspin3, unpack=True, usecols=(0,1,2,3))
        
        #plt.plot(t_ah1-t_merge, irrmass_ah1-irrmass_ah1[0], c=color[i], ls=ls[i],label=labelname)
        #plt.plot(t_ah3-t_merge, irrmass_ah3-irrmass_ah1[0], c=color[i], ls=ls[i],label=labelname)
        if 'Q2' in labelname:
            ax1.plot(t_ah1[:-100]-t_merge, irrmass_ah1[:-100]-irrmass_ah1[0], c=color[i], ls=ls[i],label=labelname)
            ax1.plot(t_ah3[530:]-t_merge, irrmass_ah3[530:]-irrmass_ah1[0], c=color[i], ls=ls[i])#,label=labelname)
        else:
            ax1.plot(t_ah1-t_merge, irrmass_ah1-irrmass_ah1[0], c=color[i], ls=ls[i],label=labelname)
            ax1.plot(t_ah3-t_merge, irrmass_ah3-irrmass_ah1[0], c=color[i], ls=ls[i])#,label=labelname)
    
        #ax2.plot(t_ah1-t_merge, np.gradient(irrmass_ah1)/np.gradient(t_ah1), c=color[i], ls=ls[i],label=labelname)
        #ax2.plot(t_ah3-t_merge, np.gradient(irrmass_ah3)/np.gradient(t_ah3), c=color[i], ls=ls[i])#,label=labelname)
    
    
        #Set all the simulations to same final time before comparing mass and spins
        area_ah3       = area_ah3    [t_ah3 < t_comp]
        irrmass_ah3    = irrmass_ah3 [t_ah3 < t_comp]
        arealrad_ah3   = arealrad_ah3[t_ah3 < t_comp]
        t_ah3          = t_ah3       [t_ah3 < t_comp]
        
        sx3   = sx3  [t_ih3 < t_comp]
        sy3   = sy3  [t_ih3 < t_comp]
        sz3   = sz3  [t_ih3 < t_comp]
        t_ih3 = t_ih3[t_ih3 < t_comp]
        
        
        bbh_irrmass_bh3 = irrmass_ah3[-1]
        bbh_spinmag_bh3 = (sx3**2 + sy3**2  + sz3**2)**0.5
        bbh_hznmass_bh3 = (bbh_irrmass_bh3**2 + 0.25*bbh_spinmag_bh3[-1]**2/bbh_irrmass_bh3**2)**0.5
        bbh_dimspin_bh3 = bbh_spinmag_bh3[-1]/bbh_hznmass_bh3**2
        
        bbh_irrmass_dict[massratio]  = bbh_irrmass_bh3
        bbh_hznmass_dict[massratio]  = bbh_hznmass_bh3
        print("%-10s \t %g \t  %g \t\t %g \t\t %g \t %g \t %g "%(labelname,bbh_irrmass_bh1, bbh_irrmass_bh3, \
                                            bbh_irrmass_bh3 - bbh_irrmass_bh3, bbh_dimspin_bh3,\
                                            bbh_hznmass_bh3,  bbh_hznmass_bh3 - bbh_hznmass_bh3))
        
        
    elif(compare_nsbh[key]['type']=='NSBH'):
       
        datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
        t_merge = compare_nsbh[key]['t_merger_src']
        t_comp = t_merge + 300.
        
        
        #if 'Q2' in labelname:
        #    datadir = os.path.join('/localdata2/bkhamesra3/simulations/Hive/NSBH/MassRatio_2/NSBH_q2_D9_a0_C0.145_M105_QLM/Summary/Data/')

        bhdiag = os.path.join(datadir,'BH_diagnostics.ah1.gp')
        t_ah, area_ah,irrmass_ah, arealrad_ah = np.loadtxt(bhdiag, unpack=True, usecols=(1,25,26,27))
        
        ihspin = os.path.join(datadir,'ihspin_hn_0.asc')
        t_ih, sx, sy, sz = np.loadtxt(ihspin, unpack=True, usecols=(0,1,2,3))
     
        mass_growth_rate = np.gradient(irrmass_ah[::8])/np.gradient(t_ah[::8])
        #plt.plot(t_ah-t_merge, irrmass_ah-irrmass_ah[0], c=color[i], ls=ls[i],label=labelname)
        ax1.plot(t_ah-t_merge, irrmass_ah-irrmass_ah[0], c=color[i], ls=ls[i],label=labelname)
        ax2.plot((t_ah-t_merge)[::8], mass_growth_rate, c=color[i], ls=ls[i],label=labelname)
        
        
        #Set all the simulations to same final time before comparing mass and spins
        area_ah       = area_ah    [t_ah<t_comp]
        irrmass_ah    = irrmass_ah [t_ah<t_comp]
        arealrad_ah   = arealrad_ah[t_ah<t_comp]
        t_ah          = t_ah       [t_ah<t_comp]
       
        sx   = sx  [t_ih<t_comp]
        sy   = sy  [t_ih<t_comp]
        sz   = sz  [t_ih<t_comp]
        t_ih = t_ih[t_ih<t_comp]
        
        
        
        nsbh_irrmass_bh1 = irrmass_ah[0]
        nsbh_spinmag_bh1 = (sx**2 + sy**2  + sz**2)**0.5
        nsbh_hznmass_bh1 = (nsbh_irrmass_bh1**2 + 0.25*nsbh_spinmag_bh1[0]**2/nsbh_irrmass_bh1**2)**0.5
        
        nsbh_irrmass_bh3 = irrmass_ah[-1]
        nsbh_spinmag_bh3 = (sx**2 + sy**2  + sz**2)**0.5
        nsbh_hznmass_bh3 = (nsbh_irrmass_bh3**2 + 0.25*nsbh_spinmag_bh3[-1]**2/nsbh_irrmass_bh3**2)**0.5
        nsbh_dimspin_bh3 = nsbh_spinmag_bh3[-1]/nsbh_hznmass_bh3**2
            
        print("%-1s \t %g \t  %g \t\t %g \t %g \t %g \t %g"%(labelname, nsbh_irrmass_bh1,nsbh_irrmass_bh3, \
                                                  (nsbh_irrmass_bh3-bbh_irrmass_dict[massratio])/bbh_irrmass_dict[massratio], nsbh_dimspin_bh3,\
                                                  nsbh_hznmass_bh3,   (nsbh_hznmass_bh3-bbh_hznmass_dict[massratio])/bbh_hznmass_dict[massratio]))
       
        
ax1.set_xlabel(r'$T/M$')
ax1.set_ylabel(r'$\delta M_{h}$')
ax2.set_xlabel(r'$T/M$')
ax2.set_ylabel(r'$dM_{h}/dt$')

ax1.set_xlim(-200,300)
ax2.set_xlim(-50,100)
#ax1.set_ylim(-2,5)
#ax2.set_ylim(-0.5,1)

ax1.legend()
ax2.legend()
#plt.axhline(y=irrmass_ah[0], ls='--', c='k')
#plt.ylim(0.65, 0.92)

ax1.grid(True)
ax2.grid(True)

plt.tight_layout()
#plt.savefig(os.path.join(figdir, 'BH_MassPlots.png'), dpi=100)
#plt.show()
plt.close()

Neutron Star Masses

The mass of neutron star rest mass remains roughly constant till the disruption at which point it drastically decreases. After this point, since the matter is consumed by the black hole, the code cannot find the surface reliably and probably leads to unreliable results.

In [ ]:
#NS Mass Plots - 
color = ['r', 'b', 'k','m','coral','forestgreen']
ls = ['-', '--', ':','-.','--',':']


plt.figure(figsize=(15,8))
for i, key in enumerate(sorted(compare_nsbh.keys())):
    
    labelname = compare_nsbh[key]['label']
    if(compare_nsbh[key]['type']=='NSBH'):
       
        datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
        t_merge = compare_nsbh[key]['t_merger_src']
        volint = os.path.join(datadir,'volume_integrals-GRMHD.asc')
        if not(os.path.exists(volint)): continue
        t, restmass = np.loadtxt(volint, unpack=True, usecols=(0,4))
        
        if compare_nsbh[key]['symmetry']:
            restmass= 2*restmass #accounting for symmetries
            
        print("%s:Initial Mass = %g, Final mass = %g"%(labelname[:4], restmass[0],restmass[-1] ))
        plt.plot(t-t_merge, restmass, c=color[i], ls=ls[i],label=labelname)
        
plt.xlabel(r'$(T - R_{ex})/M$')
plt.ylabel('Mass')

#plt.ylim(0., 2)
plt.grid(True)
plt.legend()
plt.show()
plt.close()

#NS Mass Plots - 
color = ['r', 'b', 'k','m','coral','forestgreen']
ls = ['-', '--', ':','-.','--',':']


plt.figure(figsize=(15,8))
k=0
for i, key in enumerate(sorted(compare_nsbh.keys())):
    
    labelname = compare_nsbh[key]['label']
    if(compare_nsbh[key]['type']=='NSBH'):
        k+=1
        datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
        t_merge = compare_nsbh[key]['t_merger_src']
        volint = os.path.join(datadir,'volume_integrals-GRMHD.asc')
        if not(os.path.exists(volint)): continue
        t, restmass = np.loadtxt(volint, unpack=True, usecols=(0,4))
        
        if compare_nsbh[key]['symmetry']:
            restmass= 2*restmass #accounting for symmetries
            
        print("%s:Initial Mass = %g, Final mass = %g"%(labelname[:4], restmass[0],restmass[-1] ))
        plt.subplot(2,2,k)
        plt.plot(t-t_merge, 100*(restmass[0] - restmass)/restmass[0], c=color[i], ls=ls[i],label=labelname)
        
        plt.xlabel(r'$(T - R_{ex})/M$')
        plt.ylabel(r'$\delta Mass/M_0(\%)$')

        plt.ylim(0., 2)
        plt.xlim(-t_merge, 50)
        plt.grid(True)
        plt.legend()
plt.show()
plt.close()
In [ ]:
#Accretion disk and dynamical Ejecta Rest Mass- 
color = ['r', 'b', 'k','m','coral','forestgreen']
ls = ['-', '--', ':','-.','--',':']

bbh_ah3_mass = 0


#subplt_num=1
print("System \t\t\t M(R>R_{AH})/M_* \t Accretion Disk Mass \t M_disk/M_* \t Dynamical Ejecta Mass/M_* \n ")
        
for i, key in enumerate(sorted(compare_nsbh.keys())):
    
    labelname = compare_nsbh[key]['label']
    t_merge = compare_nsbh[key]['t_merger_src']
    if (compare_nsbh[key]['type']=='BBH'):
        continue
    
    elif(compare_nsbh[key]['type']=='NSBH'):
       
        datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
   
        vi_grmhd = os.path.join(datadir,'volume_integrals-GRMHD.asc')
        time, m_ns, m_bh3, m_0d7, m_1d4, m_20, m_35  = np.loadtxt(vi_grmhd, unpack=True, usecols=(0,5,6,7,8,9,10))
        
        if compare_nsbh[key]['symmetry']:
            m_ns  = 2.*m_ns #accounting for symmetries
            m_bh3 = 2.*m_bh3 #accounting for symmetries
            m_0d7 = 2.*m_0d7 #accounting for symmetries
            m_1d4 = 2.*m_1d4 #accounting for symmetries
            m_20  = 2.*m_20 #accounting for symmetries
            m_35  = 2.*m_35 #accounting for symmetries
            
        
        if 'Q5' in labelname:tcap = 180
        else:tcap=500
        time = time-t_merge
        
        idx_cap500 = np.where(time<=tcap)
        m_ns  = m_ns [idx_cap500]
        m_bh3 = m_bh3 [idx_cap500]
        m_0d7 = m_0d7 [idx_cap500]
        m_1d4 = m_1d4 [idx_cap500]
        m_20  = m_20  [idx_cap500]
        m_35  = m_35  [idx_cap500]
        time  = time  [idx_cap500]
            
        print("%-15s \t %g \t\t %g \t\t %g\t \t  %g "%(labelname, m_1d4[-1]/m_ns[0], m_1d4[-1]-m_35[-1], \
                                               m_1d4[-1]/m_ns[0]-m_35[-1]/m_ns[0], m_35[-1]/m_ns[0] ))
        plt.figure(figsize=(15,6))
        plt.subplot(1, 2, 1)
        
        #plt.plot(time, m_0d7, c='C0', ls=':', label=labelname+'R=0.7')
        plt.plot(time, m_1d4/m_ns[0], c='C1', ls='-', label=labelname+' (R=1.5)')
        plt.plot(time, m_35/m_ns[0],  c='C2', ls='-.', label=labelname+' (R=35)')
        plt.axvline(x=t_merge, ls='-', c='k', label='Time of merger')
        plt.xlabel('Time')
        plt.ylabel(r'$M_{r>R}/M_*$')
        plt.grid(True)

        #plt.ylim(0.0, 0.3)
        #plt.xlim(t_merge-100,1100) #time[-1]+10)
        plt.legend()
        
        plt.subplot(122)
        #plt.semilogy(time, m_0d7, c='C0', ls=':', label=labelname+'R=0.7')
        plt.semilogy(time, m_1d4/m_ns[0]-m_35[-1]/m_ns[0], c='C1', ls='-', label=labelname+' (R=1.5)')
        plt.semilogy(time, m_35/m_ns[0],  c='C2', ls='-.', label=labelname+' (R=35)')
        
        plt.axvline(x=t_merge, ls='-', c='k', label='Time of merger')
        
        plt.xlabel('Time')
        plt.ylabel(r'$M_{r>r_{R}}/M_*$')
        #plt.ylim(0.9, 0.94)
        #plt.xlim(t_merge-100,1100)
        
        plt.legend()
        plt.grid(True)
        plt.show()
        plt.close()
        
    
#Accretion disk and dynamical Ejecta Rest Mass- 
color = ['r', 'b', 'k','m','coral','forestgreen']
ls = ['-', '--', ':','-.','--',':']

bbh_ah3_mass = 0
In [ ]:
#Total Rest Mass outside BH Plots - 
color = ['r', 'b', 'k','m','coral','forestgreen']
ls = ['-', '--', ':','-.','--',':']

bbh_ah3_mass = 0

#print("System \t\t\t Initial BH Mass \t Final BH Mass \t NSBH - BBH \n ")


#subplt_num=1
for i, key in enumerate(sorted(compare_nsbh.keys())):
    #print("\n\n")
    labelname = compare_nsbh[key]['label']
    t_merge = compare_nsbh[key]['t_merger_src']
    if (compare_nsbh[key]['type']=='BBH'):
        continue
    
    elif(compare_nsbh[key]['type']=='NSBH'):
       
        datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
       
        vi_grmhd = os.path.join(datadir,'volume_integrals-GRMHD.asc')
        time, m_ns = np.loadtxt(vi_grmhd, unpack=True, usecols=(0,5))
        if compare_nsbh[key]['symmetry']:
            m_ns  = 2.*m_ns #accounting for symmetries
        
        totalrestmass = os.path.join(datadir,'TotalRestMass.asc')
        if 'Q2' in labelname: totalrestmass = '/localdata2/bkhamesra3/simulations/Hive/NSBH/MassRatio_2/NSBH_q2_D9_a0_C0.145_M105_QLM/Summary/Data/TotalRestMass.asc'
        data =  np.loadtxt(totalrestmass, unpack=True, usecols=(0,1,2,4,6,8,10))
        it = data[0]
        
        print(it[2]-it[1])
        if it[2]-it[1]<30:
            
            data = (data.T[::128]).T
            it=data[0]
            print(it[2]-it[1])
            
        it, time, m_0d7, m_1d4, m_15, m_25,m_35  = data

        #print("%-15s \t %g \t \t %g \t %g"%(labelname, irrmass_ah[0],irrmass_ah[-1], irrmass_ah[-1]-bbh_ah3_mass ))
        #plt.subplot(2, 1, subplt_num)
        plt.figure(figsize=(15,8))
        #plt.plot(time, m_0d7, c='C0', ls='solid',label=labelname+'(R=0.7)')
        #plt.plot(time, m_1d4,  c='C1', ls='--',label=labelname+'(R=1.4)')
        #plt.plot(time, (m_6-m_1d4)/m_ns[0],  c='C2', ls=':', label='M(R=6) - M(R=1.4)')
        plt.plot(time, (m_25-m_1d4), c='C3', ls='-', label='M(R=25)- M(R=1.4)')
        plt.plot(time, (m_35-m_25), c='C4', ls='-.', label='M(R=35)- M(R=25)')
        
        plt.axvline(x=t_merge, ls='-', c='k', label='Time of merger')
        #plt.yscale('log')
        plt.xlabel('Time')
        plt.ylabel(r'Rest Mass /$M_*$')
        plt.title(labelname, fontsize=16)
        #subplt_num+=1
        plt.ylim(1e-6,1.2)
        plt.xlim(0,1100)
        plt.legend()
        plt.show()
        plt.close()

Rough Work

In [ ]:
#BBH Old and New comparisons

ylm1 = '/localdata2/bkhamesra3/simulations/BBH/Test/D9_q3.0_a0.0_m160/Ylm_WEYLSCAL4::Psi4r_l2_m2_r100.00.asc'
t, re, im = np.loadtxt(ylm1, unpack=True, usecols=(0,1,2))
amp = (re**2 + im**2)**0.5
plt.figure(figsize=(15,6))
plt.plot(t, re, c='C0', label='Old')
plt.plot(t, amp, c='C0',ls='--')

ylm1 = '/localdata2/bkhamesra3/simulations/BBH/Hive/BBH_q3_D9_a0_M250/BBH_q3_D9_a0_M250-all/Ylm_WEYLSCAL4::Psi4r_l2_m2_r100.00.asc'

t, re, im = np.loadtxt(ylm1, unpack=True, usecols=(0,1,2))
amp = (re**2 + im**2)**0.5
plt.plot(t, re, c='C1', label='New')
plt.plot(t, amp, c='C1',ls='--')
plt.xlim(700,950)
plt.legend()
plt.show()
plt.close()
In [ ]:
os.chdir('/localdata/bkhamesra3/research_localdata/UsefulScripts/Python_Notebooks/BowenID_Tests/NSBH/Simulation_FinalResults/MassRatio_3/')
In [ ]:
! jupyter nbconvert --to html NSBH_q3_D9_a0_C0.145_Comparison.ipynb
In [ ]:
! mv NSBH_q3_D9_a0_C0.145_Comparison.html /localdata2/bkhamesra3/website/github_pages/bkhamesra3.github.io/research/nsbh/Separation_9M/Comparisons/MassRatio3/
In [ ]:
a = np.array((1,2,3))
print((np.where(a>3) and a[-1]))
In [ ]:
# Q5 GW frequency oscillations
plt.figure(figsize=(15,8))
for i, key in enumerate(sorted(compare_nsbh.keys())):
    
    labelname = compare_nsbh[key]['label']
    if (compare_nsbh[key]['type']=='NSBH'):
        datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
        tmerge_nsbh = compare_nsbh[key]['t_merger_src']
    elif(compare_nsbh[key]['type']=='BBH'):
        datadir = compare_nsbh[key]['dir']
        continue
    
    filepath = os.path.join(datadir, 'Ylm_WEYLSCAL4::Psi4r_l2_m2_r70.00.asc')
    if not os.path.exists(filepath):
        filepath = os.path.join(datadir, 'mp_Psi4r_l2_m2_r70.00.asc')
    
    q = float((labelname.split('Q')[1]).split('-')[0])
    mtotal = 1.35*(1+q)
   
    if 'Q5' not in labelname:continue
    t22, psi4_re22, psi4_im22 = np.loadtxt(filepath, unpack=True, usecols=(0,1,2))

    #Compute GW frequency
    phi = np.unwrap(np.arctan2(psi4_im22,psi4_re22))
    amp = np.sqrt(psi4_re22**2 + psi4_im22**2)
    max_amp_idx = np.argmax(amp==np.amax(amp))
    gw_freq = np.gradient(phi)/np.gradient(t22)/2./np.pi
    
    
    #Limit the data to a small time interval during inspiral say 200-300M
    idx_fft = np.intersect1d(np.where(t22>320), np.where(t22<420))
    time_fft = t22[idx_fft]
    
    gw_freq_fit = np.poly1d(np.polyfit(t22[idx_fft], gw_freq[idx_fft], deg=2))
    tfreq_fit = np.arange(300, 500, t22[1]-t22[0])
    
    osc_fft = np.fft.fft((gw_freq[idx_fft]-gw_freq_fit(time_fft))/gw_freq_fit(time_fft))
    fft_freq = np.fft.fftfreq(np.size(time_fft), d=time_fft[2]-time_fft[1])

    plt.subplot(121)
    #plt.plot(t22, (gw_freq-gw_freq_fit(t22))/gw_freq_fit(t22),ls=ls[i], c=color[i], label=labelname)
    plt.plot(t22, gw_freq,ls='-', c=color[i], label=labelname)
    plt.plot(tfreq_fit, gw_freq_fit(tfreq_fit),ls='--', c='k', label=labelname)
    plt.xlabel(r'$Time - T_{max}$')
    plt.ylabel('Psi4 Frequency Oscillations')
    plt.xlim(200,700)
    plt.ylim(-0.02,-0.003)
    #plt.ylim(-0.1,0.1)
    
    plt.subplot(122)
    plt.plot(fft_freq/(mtotal/s1_in_Msun), osc_fft.real,ls=ls[i], c=color[i], label=labelname)
    plt.xlabel(r'$FFT Freq$')
    plt.ylabel('FFT(Psi4 Frequency)')
    plt.axvline(x=0.028)
    plt.xlim(0,2000)
    #plt.ylim(-3,2.5)
    
plt.legend()
plt.show()
plt.close()
In [ ]:
import h5py

h5file = '/localdata2/bkhamesra3/simulations/InjectionStudies/NSBH/Completed/NonSpinning/GT0014_NSBH.h5'

wf_data = h5py.File(h5file, 'r')
#$for atr_keys in wf_data.attrs.keys():
flow = wf_data.attrs['f_lower_at_1MSUN']
print(flow/8.1)
wf_data.close()